




function [U,x,t]=FDMDParabsol(f,a,b,u0,T=1,al=@(t) 0,be=@(t) 0,c=0,N=100,M=20)
%we want to solve u_t+(-u''+cu)=f in (0,T]x (a,b) with u(a,b)=al u(b,t)=be (al be constants)
%initial value u(0,x)=u0(x)
%by a FDM method on uniform mesh (x(k)=a+(k-1)*h) k=1,N+1
%u(1)=al bnd condition at a
% d/dt u(t,k) + (1/(h*h)(-u(t,k-1)+2u(t,k)-u(k+1)) +c*u(t,k)=f(t,x(k)) k=2,..,N  
%u(N+1)=be
%u(x(k),0)
% d/dt U=-AU + F + [al(t);0;0;...;0;be(t)]/(h*h)
%A= (h^(-2)d*iscrete FDM laplacian matrix + c*ID) 
%%EXAMPLE
%d/dt u-u''=0 u(0)=u(pi)=0 u(x,0)=sin(x)  sol=exp(-t)sin(x) 
%a=0;b=pi;c=0;T=1;al=be=@(t) 0;u0=@sin;f=@(t,x) zeros(size(x));N=M=20; [U,x,t]=FDMDParabsol(f,a,b,u0,T,al,be,c,N,M)
%EXAMPLE 2
%d/dt u-u''=0 u(0)=0 u(b)=exp(-t)*sin(b) u(x,0)=sin(x)  sol=exp(-t)*sin(x) 
%a=0;b=1;c=0;T=1;al=@(t) 0;be=@(t) exp(-t)*sin(b);u0=@sin;f=@(t,x) zeros(size(x));N=M=20; [U,x,t]=FDMDParabsol(f,a,b,u0,T,al,be,c,N,M)
 global ff h A x all bee
 all=al;
 bee=be;
 ff=f; %RHS function of PDE
 h=(b-a)/N;
 x=linspace(a,b,N+1);x=x';
 d=ones(N-1,1);
 V=[2*ones(N-1,1),-ones(N-1,2)];
 A=spdiags(V,[0,1,-1],N-1,N-1);
% A=sparse(2*diag(d)-diag(d(1:N-2),1)-diag(d(1:N-2),-1));
 %the last line is replaced with spdiags
 A=(1/(h*h))*A;
 A+=c*speye(size(A));
 u0x=u0(x(2:N));
 t=linspace(0,T,M+1)';
 U=lsode(@vf,u0x,t);
 U=[al(t).*ones(rows(U),1),U,be(t).*ones(rows(U),1)];
endfunction


function y=vf(u,t)
%function for lsode  
global ff h A x all bee
 N=length(x)-1;
 F=zeros(N-1,1);
 F=ff(t,x(2:N));
 %testing% [norm(F);all;bee]
 y= -A*u+F;
 y(1)+=all(t)/(h*h);
 y(end)+=bee(t)/(h*h);
endfunction
