format long
N=100;
sigma=0.2;
T=1;
r=0.05;
K=100;
alpha=4;
a=0;
b=K*alpha;
L=b-a;
h=L/(N-1);
x=0:h:L;
A=diag((2*h/3)*ones(1, N)); %Construimos A
A(1,1)=h/3;
A(N,N)=h/3;
A(1,2)=h/6;
A(N,N-1)=h/6;
for i=2:(N-1)
    A(i,i+1)=h/6;
    A(i,i-1)=h/6;
end
B=zeros(N,N); %Construimos B
B(1,1)=h/3;
B(N,N)=h/3;
for i=2:(N-1)
    xmas=x(i+1);
    xmenos=x(i-1);
    xigual=x(i);
    B(i,i+1)=(1/(3*h^2))*(xigual^3-xmas^3);
    B(i,i-1)=(1/(3*h^2))*(xmenos^3-xigual^3);
    B(i,i)=(1/(3*h^2))*(xmas^3-xmenos^3);
end

C=zeros(N,N); %Construimos C
F=sigma^2-r;
C(1,1)=-h/6;
C(N,N)=(1/(6*h^2))*(2*x(N)^3+x(N-1)^3-3*x(N)^2*x(N-1));
for i=2:(N-1)
    xmas=x(i+1);
    xmenos=x(i-1);
    xigual=x(i);
    C(i,i+1)=(-1/(6*h^2))*(-3*xmas^2*xigual-2*xigual^3+2*xmas^3+3*xigual^3);
    C(i,i-1)=(1/(6*h^2))*(-3*xmenos^2*xigual+3*xigual^3+2*xmenos^3-2*xigual^3);
    C(i,i)=(1/(6*h^2))*(xmenos^3-3*xigual^2*xmenos-xmas^3+3*xmas*xigual^2);
end

C(N,N)=C(N-1,N-1);

D=sigma^2/2*B+F*C+r*A;

Ainv=inv(A);
Matriz=-Ainv*D;


N2=N;
a2=0;
b2=T;
L2=b2-a2;
h2=L2/(N2-1);
t=0:h2:L2;
SOLUCION=zeros(N,N2);
u0 = max(x - K, 0);   % vector Nx1
u0=transpose(u0);
for j=1:N2
    SOLUCION(:,j)=expm(Matriz*t(j))*u0;
end

[X, T] = meshgrid(x,t);
surf(X, T, SOLUCION.')
shading interp        % para suavizar
colormap jet          % colores
colorbar              % referencia de colores
xlabel('x')
ylabel('t')
zlabel('V(x,t)')
hold on
SOLUCIONREAL=zeros(N,N2);
for i=1:N
    for j=1:N2
        SOLUCIONREAL(i,j)=SolExacta(x(i), t(j), K, r, sigma);
    end
end

[Xm, Tm] = meshgrid(x, t);
surf(Xm, Tm, SOLUCIONREAL.')
 
surf(Xm,Tm,SOLUCION-SOLUCIONREAL)
ERROR=SOLUCION-SOLUCIONREAL;

function V = SolExacta(S, tau, K, r, sigma)
    % Solución analítica del modelo de Black-Scholes para un call europeo
    % usando la formulación con tiempo hacia delante tau = T - t.

    % Evitar problemas numéricos cuando tau = 0
    tau = max(tau, 1e-12);

    d1 = (log(S./K) + (r + 0.5*sigma^2).*tau) ./ (sigma.*sqrt(tau));
    d2 = d1 - sigma.*sqrt(tau);

    V = S .* normcdf(d1) - K .* exp(-r*tau) .* normcdf(d2);
end