Page 124 -

b=-2*t;
c=20*ones(1,N);
y=zeros(1,N);
D=zeros(1,N);
dt=(tfin-tin)/(N-1);
u=zeros(1,N);
y(1)=3/8;
D(1)=0;
D2(1)=(1/a(1))*(-b(1)*D(1)-c(1)*y(1)+u(1));
for k=2:N
y(k)=((4*a(k)/dt^2+2*b(k)/dt+c(k))^(-1))*...
(y(k-1)*(4*a(k)/dt^2+2*b(k)/dt)+D(k-1)...
*(4*a(k)/dt+b(k))+a(k)*D2(k-1)+u(k));
D(k)=(2/dt)*(y(k)-y(k-1))-D(k-1);
D2(k)=(4/dt^2)*(y(k)-y(k-1))-(4/dt)*D(k-1)-...
D2(k-1);
end

yanal=(35*t.^4-30*t.^2+3)/8;
plot(t,y,t,yanal,'--')

As you will observe upon running this program, the numerical solution and
the analytical solution agree very well.
NOTE The above ODE is that of the Legendre polynomial of order l = 4,
encountered earlier in Chapter 2, in Pb. 2.25.

2
dP dP
(1− t 2 ) 2 l − 2t l + ll ) 1 l = 0 (4.39)
( + P
dt dt
where
P()−= (−1 ) l P t( ) (4.40)
t
l l

Homework Problem

Pb. 4.39 The above algorithms assume that the source function is continu-
ous. If it is not, we may encounter problems upon applying this algorithm
over a transition region, as will be illustrated in the following problem.

© 2001 by CRC Press LLC