% Ex c11.1a ur Heath

clear all
clf
deltax=0.05;deltat=0.0012;tsteg=50;N=20;
h=1/20; n=N-1;
x=h:h:1-h;

ettor=ones(n,1);
A=spdiags([ettor,-2*ettor,ettor],[-1,0,1],n,n)/h^2;

for i=1:n
   if i<n/2
      y0(i)=2*x(i);
   else
      y0(i)=2*(1-x(i));
   end
end
y0=y0';
Y=[y0];

for i=1:tsteg
   y1=y0+deltat*A*y0;
   Y=[Y,y1];
   y0=y1;
end;

tled=[0:tsteg]*deltat;
xled=[0,x,1];
Y=[zeros(1,tsteg+1);Y;zeros(1,tsteg+1)];
[X,T]=meshgrid(xled,tled);
mesh(X,T,Y');
hold on
title('Exempel C11.1 a');
xlabel('x - led');
ylabel('t- led');
zlabel('funktionen u(t,x)');
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Ex c11.1e ur Heath

clear all
clf
global A
deltax=0.05;N=20;
h=1/20; n=N-1;
x=h:h:1-h;

ettor=ones(n,1);
A=spdiags([ettor,-2*ettor,ettor],[-1,0,1],n,n)/h^2;

for i=1:n
   if i<n/2
      y0(i)=2*x(i);
   else
      y0(i)=2*(1-x(i));
   end
end
y0=y0';
options=odeset('RelTol',5e-2,'AbsTol',1e-2);
[tled,Y]=ode45('funkc11',[0,0.06],y0,options);
%[tled,Y]=ode45('funkc11',[0 0.06],y0);

tsteg=length(tled);
xled=[0,x,1];
Y=[zeros(1,tsteg);Y';zeros(1,tsteg)];
[X,T]=meshgrid(xled,tled);
mesh(X,T,Y');
hold on
title('Exempel C11.1 e');
xlabel('x - led');
ylabel('t- led');
zlabel('funktionen u(t,x)');
   
   
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
function y1=funkc11(t,y0)
global A
y1=A*y0;   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% Ex c11.2b ur Heath

clear all
clf
global A
N=20;
h=1/N; 
n=N+1;
x=0:h:1;

ettor=ones(n,1);
A=spdiags([ettor,-2*ettor,ettor],[-1,0,1],n,n)/h^2;
A(n,n-1)=2/h^2;
A(1,2)=2/h^2;
y0=cos(pi*x)';


[tled,Y]=ode45('funkc11',[0 0.1],y0);



tsteg=length(tled);
xled=x;

[X,T]=meshgrid(xled,tled);
subplot(2,1,1);
mesh(X,T,Y);
hold on
title('Exempel C11.2-b');
xlabel('x - led');
ylabel('t- led');
zlabel('funktionen u(t,x)');
   
   
U=exp(-pi^2*T).*cos(pi*X);
felet=max(max(abs(U-Y)))
h
subplot(2,1,2);
mesh(X,T,U-Y);
title('Felet');
xlabel('x - led');
ylabel('y - led');
zlabel('felet');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Exempel c10.1b Finite difference method.

clear all; clf;hold on;
  global A b;
M=[2,4,8,16];

for N=M
   h=1/N;
   n=N-1;
   x=[0:h:1]';
   starty=h:h:1-h;
   
 
   ett=ones(n,1);
   A=spdiags([ett,(-2-3*h^2)*ett,ett],[-1 0 1],n,n)/h^2;
   b=-([h:h:1-h].^2)';
   b(n)=b(n)+1/h^2;
   
   svar=fsolve('kurvan',starty');
   plot(x,[0;svar;1]);
   
end;
grid on;
title('Exempel c10.1b');
xlabel('t - led');
ylabel('y - led');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function z=kurvan(y)
global A b;
z=A*y-10*y.^3+b;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
