% LINEAR SHOOTING ALGORITHM 11.1
%
% To approximate the solution of the boundary-value problem
%
% Y'' =  (-2/x)Y' + 2/(x^2) Y + sin(log(x))/(x^2), 1<=X<=2, Y(1)=1, Y(2)=2:
%
% OUTPUT: Approximations W(1,I) to Y(X(I)); W(2,I) to Y'(X(I))
%         for each I=0,1,...,N.
 syms('OK', 'A', 'B', 'ALPHA', 'BETA', 'N', 'FLAG', 'NAME');
 syms('OUP', 'H', 'U1', 'U2', 'V1', 'V2', 'I', 'X', 'T');
 syms('K11', 'K12', 'K21', 'K22', 'K31', 'K32', 'K41', 'K42');
 syms('U', 'V', 'W1', 'Z', 'W2', 's', 'x');
 TRUE = 1;
 FALSE = 0;
 fprintf(1,'This is the Linear Shooting Method.\n');
 
 P = inline('-2/x','x');
 Q = inline('2/(x^2)','x');
 R = inline('sin(log(x))/(x^2)','x');
 
 %left and right endpoints 
 A = 1.0;
 B = 2.0;
 
 % Y values at left and right endpoints 
 ALPHA = 1.0;
 BETA = 2.0;
 
 %positive integer for the number of subintervals
 N = 10;
 
 OUP = 1;
 
 % W(1) is the approximation to Y(x)
 % W(2) is the approximation to Y'(x)
 fprintf(OUP, 'W(1) is the approximation to Y(x)\n\n');
 fprintf(OUP, 'W(2) is the approximation to dY(x)/dx\n\n');
 
 fprintf(OUP, 'LINEAR SHOOTING METHOD\n\n');
 fprintf(OUP, '  I  X(I)        W(1,I)      W(2,I)\n');
% STEP 1
 H = (B-A)/N;
 U1 = ALPHA;
 U2 = 0;
 V1 = 0;
 V2 = 1;
 U = zeros(2,N);
 V = zeros(2,N);
% STEP 2
 for I = 1 : N 
% STEP 3
 X = A+(I-1)*H;
 T = X+0.5*H;
% STEP 4
 K11 = H*U2;
 K12 = H*(P(X)*U2+Q(X)*U1+R(X));
 K21 = H*(U2+0.5*K12);
 K22 = H*(P(T)*(U2+0.5*K12)+Q(T)*(U1+0.5*K11)+R(T));
 K31 = H*(U2+0.5*K22);
 K32 = H*(P(T)*(U2+0.5*K22)+Q(T)*(U1+0.5*K21)+R(T));
 T = X+H;
 K41 = H*(U2+K32);
 K42 = H*(P(T)*(U2+K32)+Q(T)*(U1+K31)+R(T));
 U1 = U1+(K11+2*(K21+K31)+K41)/6;
 U2 = U2+(K12+2*(K22+K32)+K42)/6;
 K11 = H*V2;
 K12 = H*(P(X)*V2+Q(X)*V1);
 T = X+0.5*H;
 K21 = H*(V2+0.5*K12);
 K22 = H*(P(T)*(V2+0.5*K12)+Q(T)*(V1+0.5*K11));
 K31 = H*(V2+0.5*K22);
 K32 = H*(P(T)*(V2+0.5*K22)+Q(T)*(V1+0.5*K21));
 T = X+H;
 K41 = H*(V2+K32);
 K42 = H*(P(T)*(V2+K32)+Q(T)*(V1+K31));
 V1 = V1+(K11+2*(K21+K31)+K41)/6;
 V2 = V2+(K12+2*(K22+K32)+K42)/6;
 U(1,I) = U1;
 U(2,I) = U2;
 V(1,I) = V1;
 V(2,I) = V2;
 end;
% STEP 5
 W1 = ALPHA;
 Z = (BETA-U(1,N))/V(1,N);
 X = A;
 I = 0;
 fprintf(OUP, '%3d %11.8f %11.8f %11.8f\n', I, X, W1, Z);
 for I = 1 : N 
 X = A+I*H;
 W1 = U(1,I)+Z*V(1,I);
 W2 = U(2,I)+Z*V(2,I);
 fprintf(OUP, '%3d %11.8f %11.8f %11.8f\n', I, X, W1, W2);
 end;
 if OUP ~= 1 
 fclose(OUP);
 fprintf(1,'Output file %s created successfully \n',NAME);
 end;
% STEP 7