1 tic;
 2 % this method is transform from Ritz method 
 3 %is used for solving two point BVP
 4 %this code was writen by HU.D.dong in February 11th 2017
 5 %MATLAB 7.0
 6 clear
 7 clc
 8 N=50;
 9 h=1/N;
10 X=0:h:1;
11 f=inline(\'pi^2/2*sin(pi/2*x)\');
12 %以下是右端向量；
13 for i=2:N
14     fun1=@(x) f(X(i-1)+h*x).*x+f(X(i)+h*x).*(1-x);
15 fi(i-1,1)=h*quad(fun1,0,1);
16 end
17 funN=@(x) f(X(N-1)+h*x).*x;
18 fi(N,1)=h*quad(funN,0,1);
19 %以下是stiff矩阵；
20 a11=1/h+pi^2*h/12;
21 aii=2*a11;
22 aij=-1/h+pi^2*h/24;
23 A=diag(aii*ones(N,1),0)+diag(aij*ones(N-1,1),1)+diag(aij*ones(N-1,1),-1);
24 A(N,N)=a11;
25 numerical_solution=A\fi;
26 numerical_solution=[0;numerical_solution];
27 %以下是真解；
28 for i=1:length(X)
29    Accurate_solution(i,1)=sin((pi*X(i))/2)/2 - cos((pi*X(i))/2)/2 + exp((pi*X(i))/2)*((exp(-(pi*X(i))/2)*cos((pi*X(i))/2))/2 + (exp(-(pi*X(i))/2)*sin((pi*X(i))/2))/2);
30 end 
31     grid on; 
32     subplot(1,2,1);
33      plot(X,numerical_solution,\'ro-\',X,Accurate_solution,\'b^:\');
34     title(\'Numerical solutions vs Accurate solutions\');
35     legend(\'Numerical_solution\',\'Accurate_solution\');
36     subplot(1,2,2);
37       plot(X,numerical_solution-Accurate_solution,\'b x\');
38     legend(\'error_solution\');
39     title(\'error\');
40     toc;