clc
Enter the value of lower boundary: 0
clear                                                     Enter the value of upper boundary: %pi/2
function [x, y]=secondOrderODE(x1, xn, y1, yn,            Enter the value of y(l): 1
n)                                                        Enter the value of y(n): 1
h=(xn - x1)/(n -1)
x=linspace(x1,xn,n)                                        "x          y"
A=zeros(n-2,n-2)
for i=1:n-3                                                0.          1.
A(i,i)=h^2 - 2                                             0.0826735   1.0792224
                                                           0.165347    1.1510683
A(i,i + 1)= 1                                              0.2480205   1.2150469
A(i+1,i)= 1                                                0.330694    1.2707207
end                                                        0.4133675   1.3177092
                                                           0.4960409   1.3556914
A(n-2,n-2)=h^2-2
                                                           0.5787144   1.3844075
B=zeros(n-2)                                               0.6613879   1.4036613
B(1)=-y1                                                   0.7440614   1.4133212
B(n-2)= -yn                                                0.8267349   1.4133212
                                                           0.9094084   1.4036613
y(1)=y1                                                    0.9920819   1.3844075
y(2:n-1)=A\B                                               1.0747554   1.3556914
y(n)=yn                                                    1.1574289   1.3177092
return;                                                    1.2401024   1.2707207
                                                           1.3227759   1.2150469
endfunction                                                1.4054493   1.1510683
//user input part of the function                          1.4881228   1.0792224
disp("Solving second order differential equation           1.5707963   1.
using finite difference method")
n=input("Enter the number of subintervals (N):
");
x1=input("Enter the value of lower boundary: ");
xn=input("Enter the value of upper boundary: ");
y1=input("Enter the value of y(l): ");
yn=input("Enter the value of y(n): ");
[x, y]=secondOrderODE(x1, xn, y1,yn,n);
disp("x      y");
disp([x',y]);
plot2d(x,y,2)
title('finite difference method', 'fontsize',5,
'color',2,'edgecol',4,'fontname',5)
xlabel('x axis', 'fontsize',3,'color', 2,'fontname', 3)
ylabel('y axis','fontsize',3,
'color',2,'fontname',5,'rotation',0)
a=gca()
a.box ="on"
a.thickness=1
a.font_size=3
a.font_style=3
a.children.children.thickness=3
xstring(.3, 1.1,"$\frac{d^2y}{dx^2}+y=0:
y(0)=1,y(pi/2)=1$" )
d=get("hdl")
d.font_size=4