thetameth.m

various matlab programs to slove various problems

function [u,x] = thetameth(I,n,u0,f,bc,nu,theta)
nx = n(1); h = (I(2)-I(1))/nx; 
x = [I(1):h:I(2)]; t = I(3);   uold = (eval(u0))';
nt = n(2); k = (I(4)-I(3))/nt; e = ones(nx+1,1);
K = spdiags([(h/(6*k)-nu*theta/h)*e, (2*h/(3*k)+2*nu*theta/h)*e, ...
      (h/(6*k)-nu*theta/h)*e],-1:1,nx+1,nx+1);
B = spdiags([(h/(6*k)+nu*(1-theta)/h)*e, (2*h/(3*k)-nu*2*(1-theta)/h)*e, ...
      (h/(6*k)+nu*(1-theta)/h)*e],-1:1,nx+1,nx+1);
K(1,1)       = 1; K(1,2)     = 0; B(1,1)       = 0; B(1,2)     = 0;
K(nx+1,nx+1) = 1; K(nx+1,nx) = 0; B(nx+1,nx+1) = 0; B(nx+1,nx) = 0;
[L,U]=lu(K);
t = I(3); x=[I(1)+h:h:I(2)-h]; fold = (eval(f))'; 
fold = h*fold;
fold = [bc(1); fold; bc(2)];
for time = I(3)+k:k:I(4)
     t = time;     fnew = (eval(f))'; 
     fnew = h*fnew;
     fnew = [bc(1); fnew; bc(2)];
     b = theta*fnew+(1-theta)*fold + B*uold;
     y = L \ b;   u = U \ y; 
     uold = u;
end
x = [I(1):h:I(2)];
return