Matlab to solve finite difference for Ode

% Finite difference method

% Approximate the solution of y"=(-2/x)y'+(2/x^2)y+ sin(lnx)/x^2

% for 1<=x<=2 with y(1)=1 and y(2)=2.

p = @(x) (-2/x);

q = @(x) (2/x^2);

r = @(x) (sin(log(x))/x^2);

aa = 1; bb = 2; alpha = 1; beta = 2; n=29;      % h = (bb-aa)/(n+1);   h=0.1

fprintf('   x           w   \n');

h = (bb-aa)/(n+1);

a = zeros(1,n+1);

b = zeros(1,n+1);

c = zeros(1,n+1);

d = zeros(1,n+1);

l = zeros(1,n+1);

u = zeros(1,n+1);

z = zeros(1,n+1);

w = zeros(1,n+1);

x = aa+h;

a(1) = 2+h^2*q(x);

b(1) = -1+0.5*h*p(x);

d(1) = -h^2*r(x)+(1+0.5*h*p(x))*alpha;

m = n-1;

for i = 2 : m

x = aa+i*h;

a(i) = 2+h^2*q(x);

b(i) = -1+0.5*h*p(x);

c(i) = -1-0.5*h*p(x);

d(i) = -h^2*r(x);

end

x = bb-h;

a(n) = 2+h^2*q(x);

c(n) = -1-0.5*h*p(x);

d(n) = -h^2*r(x)+(1-0.5*h*p(x))*beta;

l(1) = a(1);

u(1) = b(1)/a(1);

z(1) = d(1)/l(1);

for i = 2 : m

l(i) = a(i)-c(i)*u(i-1);

u(i) = b(i)/l(i);

z(i) = (d(i)-c(i)*z(i-1))/l(i);

end

l(n) = a(n)-c(n)*u(n-1);

z(n) = (d(n)-c(n)*z(n-1))/l(n);

w(n) = z(n);

for j = 1 : m

i = n-j;

w(i) = z(i)-u(i)*w(i+1);

end

i = 0;

fprintf('%5.4f    %11.8f\n', aa, alpha);

for i = 1 : n

x = aa+i*h;

fprintf('%5.4f    %11.8f\n', x, w(i));

end

i = n+1;

fprintf('%5.4f    %11.8f\n', bb, beta);am trying to run d code but it keeps writing Error: Invalid text character. The text ' ' contains an unsupported non-ASCII whitespace character. Please where am I doing wrong?

Answer 1

[Viranch Patel]

I tried running the code on my machine and it's working fine. You can paste this code to your machine and see if its working. I am attaching the output too. You might be using some wrong character for any particular character which is not supported.

% Finite difference method
% Approximate the solution of y"=(-2/x)y'+(2/x^2)y+ sin(lnx)/x^2
% for 1<=x<=2 with y(1)=1 and y(2)=2.
p = @(x) (-2/x);
q = @(x) (2/x^2);
r = @(x) (sin(log(x))/x^2);
aa = 1; bb = 2; alpha = 1; beta = 2; n=29;      % h = (bb-aa)/(n+1);   h=0.1
fprintf('   x           w   \n');
h = (bb-aa)/(n+1);
a = zeros(1,n+1);
b = zeros(1,n+1);
c = zeros(1,n+1);
d = zeros(1,n+1);
l = zeros(1,n+1);
u = zeros(1,n+1);
z = zeros(1,n+1);
w = zeros(1,n+1);
x = aa+h;
a(1) = 2+h^2*q(x);
b(1) = -1+0.5*h*p(x);
d(1) = -h^2*r(x)+(1+0.5*h*p(x))*alpha;
m = n-1;
for i = 2 : m
    x = aa+i*h;
    a(i) = 2+h^2*q(x);
    b(i) = -1+0.5*h*p(x);
    c(i) = -1-0.5*h*p(x);
    d(i) = -h^2*r(x);
end
x = bb-h;
a(n) = 2+h^2*q(x);
c(n) = -1-0.5*h*p(x);
d(n) = -h^2*r(x)+(1-0.5*h*p(x))*beta;
l(1) = a(1);
u(1) = b(1)/a(1);
z(1) = d(1)/l(1);
for i = 2 : m
    l(i) = a(i)-c(i)*u(i-1);
    u(i) = b(i)/l(i);
    z(i) = (d(i)-c(i)*z(i-1))/l(i);
end
l(n) = a(n)-c(n)*u(n-1);
z(n) = (d(n)-c(n)*z(n-1))/l(n);
w(n) = z(n);
for j = 1 : m
    i = n-j;
    w(i) = z(i)-u(i)*w(i+1);
end
i = 0;
fprintf('%5.4f    %11.8f\n', aa, alpha);
for i = 1 : n
    x = aa+i*h;
    fprintf('%5.4f    %11.8f\n', x, w(i));
end
i = n+1;
fprintf('%5.4f    %11.8f\n', bb, beta);

Output:

   x           w   
1.0000     1.00000000
1.0333     1.03067949
1.0667     1.06155254
1.1000     1.09262608
1.1333     1.12390423
1.1667     1.15538887
1.2000     1.18708018
1.2333     1.21897697
1.2667     1.25107701
1.3000     1.28337728
1.3333     1.31587416
1.3667     1.34856355
1.4000     1.38144105
1.4333     1.41450201
1.4667     1.44774163
1.5000     1.48115505
1.5333     1.51473732
1.5667     1.54848354
1.6000     1.58238883
1.6333     1.61644834
1.6667     1.65065735
1.7000     1.68501118
1.7333     1.71950528
1.7667     1.75413522
1.8000     1.78889666
1.8333     1.82378540
1.8667     1.85879736
1.9000     1.89392857
1.9333     1.92917520
1.9667     1.96453354
2.0000     2.00000000