function [u]=fixeddiffusion(h,N,dt,t0,theta,u0,uN,uinitial,x0,xN,uexact)
x=zeros(N+1,1);
x(1)=x0;
x(N+1)=xN;
uvec=zeros(N+1,1);
uvec(1)=u0;
uvec(N+1)=uN;
A=zeros(N-1,N-1);
a=-2/(h^2);
for i=1:N-1
A(i,i)=a;
end
b=(1/(h^2));
for i=1:N-2
j=i+1;
A(i,j)=b;
end
for i=2:N-1
k=i-1;
A(i,k)=b;
end
b=zeros(N-1,1);
for i=1:N-1
b(i)=uinitial(i+1);
end
uu=A\b;
for i=1:N-1
uu(i)=uvec(i+1);
end
u=zeros(N+1,1);
u(1)=u0;
for i=2:N
u(i)=u(i-1)+dt*((1-theta)*uvec(i-1)+theta*uvec(i));
end
t=zeros(N+1);
for i=1:N+1
t(i)=t0+(i-1)*dt;
end
ureal=zeros(N+1,1);
for i=1:N+1
ureal(i)=uexact(x(i),t(i));
end
figure(1)
plot(t,u,'.',t,ureal,'--');
legend('Approx:','True:');
error=zeros(N+1,1);
error=abs(u-ureal);
norm=zeros(N+1,1);
norm(i)=symsum((h*dt(abs(error(i)^2)))^(1/2),error(i),1,i);
end
uinitial=@(x,t)sin(pi*x);
u0=0;
uN=0;
h=0.1;
N=1/0.1;
dt=0.0005;
theta=0;
x0=0;
xN=1;
uexact=@(x,t)exp((-pi^2)*t)*sin(pi*x);
[u]=fixeddiffusion(h,N,dt,theta,u0,uN,uinitial,x0,xN,uexact)
