//Write Functions.

function [D1,D2] = matder(h,N)
  h2 = 2*h;
  hs = h^2;
  D1 = zeros(N,N);
  D2 = zeros(N,N);
  for i = 2:N-1
    D1(i,i-1)=-1;
    D1(i,i+1)=1;
    D1(1,1) = -3;
    D1(1,2) = 4;
    D1(1,3) = -1;
    D1(N,N) = 3;
    D1(N,N-1) = -4;
    D1(N,N-2) = 1;
    D2(i,i-1)=1;
    D2(i,i)=-2;
    D2(i,i+1)=1;
  end
  D1 = D1/h2;
  D2 = D2/hs;
endfunction

function t = tf(r,q)
  if q(6) == 1
    t = 1;
  else
    t = (q(1)*q(2))/r;
  end
endfunction

function C=Cf(r,q)
  if q(6) == 1
    C = 1/r;
  else
    dtr = -(q(1)*q(2))/(r^2);
    C = (1-((r/t)*dtr))/r;
  end
endfunction

function D=Df(r,q)
  if q(6) == 1
    D = -1/(r^2);
  else
    dtr = -(q(1)*q(2))/(r^2);
    D = ((((q(3)*r)/t)*dtr)-1)/(r^2);
  end
endfunction

function E = Ef(r,q)
  if q(6) == 1
    t = 1;
  else
    t = (q(1)*q(2))/r;
  end
  E = -((3+q(3))*(q(4)^2)*q(5)*r*t);
endfunction




//Input.
a = 4;
b = 12;
N = 6;
ca = 0;
da = 1;
ea = 0;
cb = 0;
db = 1;
eb = 0;
ta= 2;
mu= 0.3;
omega= ((4000*%pi)/30);
rho= 0.28;
r = linspace(a,b,N);
h = r(2)-r(1);
selector = 1
q = [ta,a,mu,omega,rho,selector];  //[ta,a,mu,omega,rho,selector]

//Data Setup.

C = zeros(N,N);
D = zeros(N,N);
E = zeros(N,1);
for j = 2:N-1
  C(j,j) = Cf(r(j),q);
  D(j,j) = Df(r(j),q);
  E(j) = Ef(r(j),q);
end

//Call the matder for D1 and D2.

[D1,D2] = matder(h,N);

//Assemble A.
A = D2 + C*D1 + D;

//Apply Boundary Conditions
A(1,1) = 1;  //E(1) = 0
A(N,N) = 1;  //E(N) = 0

//Solve for Y.

Y = A\E;

//Finding st and sr.
sr = Y/(r'*t);
D1b = D1;
D1b(1,1) = da - 3*ca/(2*h);
D1b(1,2) = 2*ca/h;
D1b(1,3) = -ca/(2*h);

D1b(N,N-2) = cb/(2*h);
D1b(N,N-1) = -2*cb/h;
D1b(N,N) = db + 3*cb/(2*h);

E(1) = ea;
E(N) = eb;

Yp = D1b * Y;
st = Yp + q(5)*(q(4)^2)*(r'^2);

// Plot
Y
A
E
sr
D1b
Yp
st
plot(r,st,'*',r,sr,'*')