clear all
a=1; b=3;
h=0.1;
x=a:h:b; %El vector x tiene N+1 elementos
%nuestro vector x es ro.
N=(b-a)/h;
xx=x(2:N); %quitamos primer y último elemento de x (DIRICHLET)
xx=xx';
ua=3;ub=0; %condiciones de contorno
U0=(exp(-8*(xx-2).^2)+3-xx);
% matriz K y B
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
B=1/(2*h)*(diag(ones(N-2,1),-1)-diag(ones(N-2,1),1));
B=diag(1/xx)*B;
% término F y valor inicial
F=0*xx;
F(1)=F(1)+ua/h^2-2*ua/h;
F(end)=F(end)+ub/h^2-2*ub/h;
% discretización en t
t0=0;tM=10;
k=h/4; % paso en t.
t=t0:k:tM;
M=length(t)-1; %número de subintervalos
sol(:,1)=U0;
for i=1:M
 sol(:,i+1)=(eye(size(K))+k/2*(K+B))\(sol(:,i)+k/2*(F-(K+B)*sol(:,i)+F)); %Metodo del trapecio
end
%incluimos condiciones DIRICHLET
UA=ua*ones(1,length(t));
UB=ub*ones(1,length(t));
sol=[UA;sol;UB];
% Gráfico
[Mt,Mx]=meshgrid(t,x);
mesh(Mx,Mt,sol)
