clear all
syms Y(z)
syms z P Q w Ra_np Nd NB 
syms TU TL Ra_s Kc KC Le
I=7;
a   = sym('a'   ,[1 I]);
b   = sym('b'   ,[1 I]);
c   = sym('c'   ,[1 I]);
dW  = sym('dW'  ,[1 I]);
dT  = sym('dT'  ,[1 I]);
dPS = sym('dPS' ,[1 I]);
SW  = sym('SW'  ,[1 I]);
ST  = sym('ST'  ,[1 I]);
SPS = sym('SPS' ,[1 I]);
format longG
clc
%--------------------------------------------------------
DY    = diff(Y,1);
D2Y   = diff(Y,2);
%--------------------------------------------------------
Case  = 1;
w1    = 0;
w2    = 10;

Kc    = 1.0;
NB    = 0.01;
Le    = 50;
P     = 0.10;
Ra_np = 0.50;
Nd    = 2.0;
Q     = 100.0;
TU    = 0.0;
TL    = 0.0;
%------------------------------------------ CASE I -----------------------:
if Case == 1
    Tb  = dsolve(diff(Y, 2)+Ra_s ==0, Y(0)==TL, Y(1)==TU);
    D1Tb  = diff(Tb, 1, z);
    D1Tb0 = subs(D1Tb,z,0);
    D1Tb1 = subs(D1Tb,z,1);
    KC=Kc;
    if Kc==0 
        PSb = dsolve(diff(Y, 2)+Nd*diff(Tb, 2)-KC*Le*Y-KC*Le==0, DY(0)+Nd*D1Tb0==0, Y(0)==0);
    else
        PSb = dsolve(diff(Y, 2)+Nd*diff(Tb, 2)-KC*Le*Y-KC*Le==0, DY(0)+Nd*D1Tb0==0, DY(1)+Nd*D1Tb1==0);
    end
    WW  = z*(1-2*z+z^2);
    TT  = (1-z);
    PPS = -Nd*(1-z);
end
if Case == 2
    Tb  = dsolve(diff(Y, 2)+Ra_s ==0, DY(0)==0, Y(1)==TU);
    D1Tb  = diff(Tb, 1, z);
    D1Tb0 = subs(D1Tb,z,0);
    D1Tb1 = subs(D1Tb,z,1);
    KC=Kc;
    if Kc==0 
        PSb = dsolve(diff(Y, 2)+Nd*diff(Tb, 2)-KC*Le*Y-KC*Le==0, DY(1)+Nd*D1Tb1==0, Y(0)==0);
    else
        PSb = dsolve(diff(Y, 2)+Nd*diff(Tb, 2)-KC*Le*Y-KC*Le==0, DY(0)+Nd*D1Tb0==0, DY(1)+Nd*D1Tb1==0);
    end
    WW  = z*(1-2*z+z^2);
end
if Case == 3
    Tb  = dsolve(diff(Y, 2)+Ra_s ==0, DY(1)==0, Y(0)==TL);
    D1Tb  = diff(Tb, 1, z);
    D1Tb0 = subs(D1Tb,z,0);
    D1Tb1 = subs(D1Tb,z,1);
    KC=Kc;
    if Kc==0 
        PSb = dsolve(diff(Y, 2)+Nd*diff(Tb, 2)-KC*Le*Y-KC*Le==0, DY(1)+Nd*D1Tb1==0, Y(0)==0);
    else
        PSb = dsolve(diff(Y, 2)+Nd*diff(Tb, 2)-KC*Le*Y-KC*Le==0, DY(0)+Nd*D1Tb0==0, DY(1)+Nd*D1Tb1==0);
    end
    WW  = z*(1-2*z+z^2);
end
PSb = collect(PSb);
D1PSb = diff(PSb, 1, z);
% plotTb;
% plotPSb;
W   = 0;
T   = 0;
PS  = 0;
if Case == 1
    for i=1:I
        W  = W  + a(i)*z^i*WW;
        T  = T  + b(i)*z^i*TT;
        PS = PS + c(i)*z^i*PPS;
        dW(i)  = z^i*WW;
        dT(i)  = z^i*TT;
        dPS(i) = z^i*PPS;
    end
elseif Case == 2
    for i=1:I
        W  = W  + a(i)*z^i*WW;
        T  = T  + b(i)*(1-z^(i+1));
        PS = PS + c(i)*(1-z^(i+1));
        dW(i)  = z^i*WW;
        dT(i)  = 1-z^(i+1);
        dPS(i) = 1-z^(i+1);
    end
elseif Case == 3
    for i=1:I
        W  = W  + a(i)*z^i*WW;
        T  = T  + b(i)*(z^(i+1)*exp(-(i+1)*z));
        PS = PS + c(i)*(z^(i+1)*exp(-(i+1)*z));
        dW(i)  = z^i*WW;
        dT(i)  = z^(i+1)*exp(-(i+1)*z);
        dPS(i) = z^(i+1)*exp(-(i+1)*z);
    end
end
D4W = diff(W, 4, z);
D2W = diff(W, 2, z);
D1W = diff(W, 1, z);
F1  = D4W-(2*w^2+Q+P)*D2W+(w^4+P*w^2)*W-w^2*(T-Ra_np*PS);
for i=1:I
    SW(i)=int(F1*dW(i),z,0,1);
end

D1T  = diff(T, 1, z);
D2T  = diff(T, 2, z);
D1PS = diff(PS, 1, z);
D2PS = diff(PS, 2, z);

F2  = D2T+NB*(D1PSb+2*Nd*D1Tb)/Le*D1T-w^2*T+(NB*D1Tb)/Le*D1PS-D1Tb*W;
for i=1:I
    ST(i)=int(F2*dT(i),z,0,1);
end

F3  = 1/Le*(D2PS-w^2*PS)+Nd/Le*(D2T-w^2*T)-D1PSb*W-Kc*PS;
for i=1:I
    SPS(i)=int(F3*dPS(i),z,0,1);
end

[A,b]=equationsToMatrix([SW(1:I)==0,ST(1:I)==0,SPS(1:I)==0],[a(1:I),b(1:I),c(1:I)]);
A=4e3*A;
minRas = 100000000;
minw   = 4.0;
gr     = (sqrt(5) - 1) / 2;
aa      = 0;
bb      = 30;
cc      = bb - gr * (bb - aa);
dd      = aa + gr * (bb - aa);
tol    = 2e-3;
fileID = fopen('OUT.txt','w');
while abs(cc - dd) > tol
    abs(cc - dd)
    Amatc       = subs(A,{w},{cc});
    Amatd       = subs(A,{w},{dd});
    detAc       = vpa(det(Amatc));
    detAd       = vpa(det(Amatd));
    sol_Ra_sc   = solve(detAc==0,Ra_s);
    sol_Ra_sd   = solve(detAd==0,Ra_s);
    sol_Ra_s_rc = sol_Ra_sc(find(imag(sol_Ra_sc)==0));
    sol_Ra_s_rd = sol_Ra_sd(find(imag(sol_Ra_sd)==0));
    sol_Ra_s_mc = sol_Ra_rc(find(sol_Ra_s_rc > 0));
    sol_Ra_s_md = sol_Ra_rd(find(sol_Ra_s_rd > 0));
    Sol_Rasc     = min(eval(sol_Ra_s_mc));
    Sol_Rasd     = min(eval(sol_Ra_s_md));
    if Sol_Rasc < Sol_Rasd
        bb = dd;
        dd = cc;
        cc = bb - gr * (bb - aa);
    else
        aa = cc;
        cc = dd;
        dd = aa + gr * (bb - aa);
    end
    fprintf(fileID,'%10.5f %10.5f %10.5f\r\n',aa,bb,Sol_Rasc);
end
minRas = Sol_Rasc
minw   = aa
fclose(fileID);