clear all
close all
% du/dt = d^2u/dx^2+5exp(t)sin(x)cos(x)
% u(0,x)=sin(x)cos(x)
% u(t,0)=u(t,pi)=0
% soluzione u(t,x)=exp(t)sin(x)cos(x)
theta = 1/2;
% discretizzazione spaziale
tstar  = 1;
mrange = [10:10:50];
timesteps = 2000;
k = tstar/timesteps;
counter = 0;
for m = mrange
  counter = counter+1;
  x = linspace(0,pi,m)';
  u0 = sin(x).*cos(x);
  h = pi/(m-1);
  A = toeplitz(sparse([-2/h^2,1/h^2,zeros(1,m-2)]));
  A(1,:) = 0;
  A(m,:) = 0;
  b = sin(x).*cos(x);
  b(1) = 0;
  b(m) = 0;
  I = speye(m);
  u = u0;
  t = 0;
  [L,U,P] = lu(I-(1-theta)*k*A);
  for i = 1:timesteps
    u = U\(L\(P*((I+theta*k*A)*u+k*theta*5*exp(t)*b+k*(1-theta)*5*exp(t+k)*b)));
    t = t+k;
  end  
  errors(counter) = norm(u-exp(tstar)*u0,inf);
end
loglog(mrange,errors,'*',mrange,(mrange/mrange(1)).^(-2)*errors(1))
xlabel('m')
ylabel('error')
legend('spatial error','order 2')
% discretizzazione temporale
m = 2000;
x = linspace(0,pi,m)';
u0 = sin(x).*cos(x);
h = pi/(m-1);
A = toeplitz(sparse([-2/h^2,1/h^2,zeros(1,m-2)]));
A(1,:) = 0;
A(m,:) = 0;
I = speye(m);
b = sin(x).*cos(x);
b(1) = 0;
b(m) = 0;
timestepsrange = [10:5:25];
counter = 0;
for timesteps = timestepsrange
  counter = counter+1;
  k = tstar/timesteps;
  u = u0;
  t = 0;
  [L,U,P] = lu(I-(1-theta)*k*A);
  for i = 1:timesteps
    u = U\(L\(P*((I+theta*k*A)*u+k*theta*5*exp(t)*b+k*(1-theta)*5*exp(t+k)*b)));
    t = t+k;
  end  
  errort(counter) = norm(u-exp(tstar)*u0,inf);
end  
figure
loglog(timestepsrange,errort,'*',timestepsrange,...
(timestepsrange/timestepsrange(1)).^(-2*theta-1)*errort(1))
legend('temporal error',sprintf('order %d',1+2*theta))
xlabel('timesteps')
ylabel('error')