
%FDM solver for u_t-u_{xx}=f in (a,b) with Dirichlet bc and u(0,x)=u0(x)
%ODE solver - lsode()


function [U,x,t]=FDM1DParab17(a,b,T,N,M,u0,uL,uP,f)
%a,b - ends of interval 
%T final time
%N,M no of space/time mesh points
%u0 a function handle to intial value
%uL,uP- f. handles for time functions with left/right Ditichlet bnd
%f - f. handle to the rhs of parabolic equation
%EXAMPLE
%uL=uP=@(t) zeros(size(t)); u0=@sin;a=0;b=pi; f=@(x) zeros(size(x)); N=20;M=10; T=1;
%U=FDM1DParab17(a,b,T,N,M,u0,uL,uP,f);plot(U')
global A F h UL UP x 
U=zeros(M+1,N);
x=linspace(a,b,N+1)';
h=(b-a)/N;
t=linspace(0,T,M+1);
A=sparse(2*diag(ones(N-1,1)) -diag(ones(N-2,1),-1)-diag(ones(N-2,1),1));
F=f;UL=uL;UP=uP;


xi=x(2:N);

U=lsode(@rhs,u0(xi),t);
Ua=uL(t)';Ub=uP(t)';
U=[Ua,U,Ub];
end


function y=rhs(U,t)
%vector field for parabolic problem
global A F h UL UP x 
 NN=length(x)-1;
 FF=F(t,x(2:NN));
 FF(1)=FF(1)+ UL(t)/(h*h);
 FF(end)=FF(end)+UP(t)/(h*h);
 y=-(1/(h*h))*(A*U)+FF;

end