function [x,t,Z]=CrNich(u0,ga,gb,ab,nx,nt,dt)
% [x,t,Z]=CrNich(u0,ga,gb,ab,nx,nt,dt)
%  LämpöYhtälö Crank-Nicholsonin (implisiittisellä) menetelmällä
%  23.3.05 
% u0 - funktio x->alkuarvo
% ga - funktio t->reunaehto kun x=a
% gb - funktio t->reunaehto kun x=b
% ab - 2-vektori, x kuuluu välille [a,b]
% nx - x hilapisteiden lkm
% nt - aika-askelten lkm
% dt - aika-askeleen pituus
%
a=ab(1); b=ab(2);
t=0:dt:nt*dt;
x=linspace(a,b,nx+2);
xs=x(2:end-1)';        % Sisäsolmut.
dx=(b-a)/(nx+1);
r=dt/dx^2;
e=ones(nx,1);
A=spdiags([-e 2*e -e],-1:1,nx,nx);
M=2*speye(nx,nx)+r*A; 
[Ala,Yla]=lu(M);  % Talletetaan vain ala-ja yläkolmiomatriisit yhtälösysteemin
                  % ratkaisua varten.
M=2*speye(nx,nx)-r*A; % Uusi M, jolla kerrotaan oikean puolen vm
U=zeros(nx,nt+1);
U(:,1)=u0(xs);
for m=1:nt
    uk=M*U(:,m); % 'lastiksi' M*edellinen 
    uk(1)=uk(1)+r*(ga(t(m))+ga(t(m+1)));   %reunaehdot 'lastiin'
    uk(end)=uk(end)+r*(gb(t(m))+gb(t(m+1)));   
    z=Ala\uk; 
    U(:,m+1)=Yla\z;
end
Z=[ga(t); U ;gb(t)]'; % Liitetään reunat ja transponoidaan.
