%   Author: Reija Ristilä <reija@hut.fi>
%   Last modified: 1999-04-01 16:31:04 EEST

% Problem:
% -u''(x)=f(x), u(0)=u(1)=0, f(x)=(3x+x^2)exp(x) 

m=13;
x=linspace(0,1,m);
h=x(2:m)-x(1:m-1)
inind=2:m-1;
xin=x(inind);%%hin=h(inind);
d=1./h(1:m-2)+1./h(2:m-1);
ud=-1./h(2:m-2);ld=-1./h(1:m-3);
A=spdiags([[ud,0]', d', [0,ld]'],-1:1,m-2,m-2);
full(ans)

% function y=ff(x)
% y=(3*x+x.^2).*exp(x);

% function y=xxf(x)
% y=x.*ff(x);

clear b
for i=1:m-2;
  b(i)=(quad8('xxf',x(i),x(i+1))-x(i)*quad8('ff',x(i),x(i+1)))/h(i);
  b(i)=b(i)+(x(i+2)*quad8('ff',x(i+1),x(i+2)) - quad8('xxf',x(i+1),x(i+2)))/h(i+1);
end;

% numeerinen
v=A\b';v0=[0;v;0];
% eksakti   
u=x.*(1-x).*exp(x);  

subplot(1,2,1),plot(x,u,x,v0,'--o'),title('f(x)=(3x+x^2) exp(x)') 
subplot(1,2,2),plot(x,u'-v0,'o-'),title('u-v')

% -------------------------------------------

% Problem:
% -u''(x)=f(x), u(0)=u(1)=0, f(x)=10*(1-10*(x-1/2)^2)*exp(-5*(x-1/2)^2);

%function y=ff2(x)
%y=10*(1-10*(x-1/2).^2).*exp(-5*(x-1/2).^2);

%function y=xxf2(x)
%y=x.*ff2(x);

clear b
for i=1:m-2;
  b(i)=(quad8('xxf2',x(i),x(i+1))-x(i)*quad8('ff2',x(i),x(i+1)))/h(i);
  b(i)=b(i)+(x(i+2)*quad8('ff2',x(i+1),x(i+2)) - quad8('xxf2',x(i+1),x(i+2)))/h(i+1);
end;

% numeerinen
v=A\b';v0=[0;v;0];
% eksakti
u=exp(-5*x.^2+5*x-5/4)-exp(-5/4);

subplot(1,2,1),plot(x,u,x,v0,'--o'),title('Harj5, prob 1') 
subplot(1,2,2),plot(x,u'-v0,'o-'),title('u-v')

