I want to use fzero to find the root of a transcendental function r. I have problem on line 33-37 and 52-63

%% initialization
clear 
clc
close all

%% Model parameters
global k1 eta1 alpha3 gamma1 d N Phi xi h C q qstar Theta OddAsymptotes m ...
    qq r1 r2 alpha aa b Efun F B Ainv r v theta_sol Vfun v_sol p ModifiedOddAsym ...
    qstarTemporary theta_init z
k1 = 6*10^(-12);            %  elastic constant
eta1 = 0.0240;              %  viscosity
xi = 0.585;                  %  activity parameter
alpha3 = -0.001104;         %  viscosity
gamma1 = 0.1093;            %  viscosity
Theta = 0.0001;
theta_init = Theta.*sin(pi*z/d);
d = 0.0002;
N = 20; 
% Simplified parameters
alpha=1-alpha3^2/gamma1/eta1;
C = (alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);
% Plotting r(q)
q=0:0.01:(5*pi);
r=q-(1-alpha).*tan(q)+ C.*tan(q);

Error using .*
Arrays have incompatible sizes for this operation.

figure
plot(q,r)
axis([0 5*pi -20 20])
% Assign g for q and plot r, and compute the special asymptotes
syms C
for q=0:N
    qstar(i) = solve(1/C == 0);
end
% Compute Odd asymptote
for j=1:N
    OddAsymptotes(j)= 1/2*(2*j + 1)*pi;
end
ModifiedOddAsym = [abs(OddAsymptotes) abs(qstar)];
qstarTemporary = min(ModifiedOddAsym);
% Compute the rest of the Roots
for j=1:N
    OddAsymptotes(j)= 1/2.*(2*j + 1)*pi;
end
% Compute all symptotes
AllAsymptotes = [sort(OddAsymptotes) sort(qstar)]; 
q0 =0;
for q = 0.001:AllAsymptotes(1);
    r1 = fzero(r, q0)
end
for i=1:N
    for q =AllAsymptotes(i):AllAsymptotes(i+1);
        r2 = fzero(r, q0)
    end
end 

% Compute all roots
gg = [r1 r2];
% Remove the repeated roots if any & redefine n 
qq = unique(gg);
m = length(qq);
% Compute the time constants tan_n
for i=1:m
    p(i) = (gamma1*alpha)./(4*k1*qq(i)^2/d^2-alpha3*xi/eta1);
end 
% Compute integral coefficients for off-diagonal elements θ_n matrix
for i=1:m 
    for j= i+1:m
        for z=0:d
            b(i, j) = int(Efun(i).*Efun(j));
            b(j, i) = b(i, j); 
            aa(i, j) = b(i, j);
            aa(j, i)= b(j, i);
        end 
    end
end
% Calculate integral coefficients for diagonal elements in theta_n matrix
for i=1:m
    for z=0:d
        aa(i, i)= int(Efun(i).^2);
    end 
end
% Calculate integrals for RHS vector
for i=1:m 
    for z=0:d
        F(i) = int(theta_init.*Efun(i));
    end 
end 
% Create matrix A and RHS vector B
for i=1:m
    for j= 1:m
        aa(i, j)
    end
end
for i=1:m
    B = F(i)
end
% Calculate inverse of A and solve A*x=B
Ainv = inv(A);
r = B*Ainv;
% Define Theta(z,t)
for i=1:m
    theta_sol(i) = sum(r(i).*Efun(i).*exp(-t/p(i)));
end
% Compute v_n for n times constant

for i=1:m 
    v(i) = (-2*k1*alpha3*qq(i))*(1/(d^2*eta1*alpha*gamma1))+ alpha3^2*xi/(2*(eta1)^2*qq(i)*alpha*gamma1)+xi/(2*eta1*qq(i));
end
% Compute v(z,t) from initial conditions
for i=1:m
    Vfun(i) = d*sin(qq(i)-2*qq(i)*z/d)+(2*z-d)*sin(qq(i));
end

% Define v(z,t)
for i=1:m
    v_sol = sum(v(i).*r(i)*Vfun(i).*exp(-t/p(i)));
end 

[Steven Lord]

What does "have problem" mean in this context?

• Do you receive warning and/or error messages? If so the full and exact text of those messages (all the text displayed in orange and/or red in the Command Window) may be useful in determining what's going on and how to avoid the warning and/or error.
• Does it do something different than what you expected? If so, what did it do and what did you expect it to do?
• Did MATLAB crash? If so please send the crash log file (with a description of what you were running or doing in MATLAB when the crash occured) to Technical Support so we can investigate.

Also please indicate specifically which lines are lines 33-37 and 52-63 in the code you posted above.

[University Glasgow]

I receive the following erros:

Array indices must be positive integers or logical values. Error in Case1_Ijuptilk (line 36)

qstar(i) = solve(1/C == 0);

Error using fzero :FUN must be a function, a valid character vector expression, or an inline function object. Error in Case1_Ijuptilk (line 56)

r1 = fzero(r, q0);

Error using fzero :FUN must be a function, a valid character vector expression, or an inline function object. Error in Case1_Ijuptilk (line 61)

r1 = fzero(r, q0);

[Torsten]

Array indices must be positive integers or logical values. Error in Case1_Ijuptilk (line 36)

qstar(i) = solve(1/C == 0);

i must be q. But what variable do you want to solve for ? C is quadratic in q - so it permits at most two solutions.

r1 = fzero(r, q0);

r in your code is an array of values for q =0:0.01:(5*pi). "fzero" needs a function as first argument.

Define

alpha=1-alpha3^2/gamma1/eta1;
C = @(q) (alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);
r = @(q) q-(1-alpha).*tan(q)+ C(q).*tan(q);

[University Glasgow]

I want to solve for 1/C for q from 0 to N. I tried something like ths

C = @(q) (alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);

for q=1:N

qstar(q) = solve(1/C);

end

[Torsten]

I don't understand what you mean by "I want to solve for 1/C". Do you mean: "I want to evaluate 1/C for q = 1,...,N" ?

[University Glasgow]

Thank you for your help. Please hold on for now. I'm new to matlab.

Can help me to check why the code below print same values of q

d = 0.2;
N =10;
fun = @(x, d) d.*sin((i)*x);
for i=1:N
    q = integral(@(x) fun(x, d) ,0,1,'ArrayValued',true)
end

[Torsten]

I do not even know what the code should give as result.

Do you mean

d = 0.2;
N = 10;
fun = @(x, i, d) d.*sin(i*x);
for i=1:N
    q(i) = integral(@(x) fun(x, i, d) ,0,1);
end
q

q = 1×10

0.0919 0.1416 0.1327 0.0827 0.0287 0.0013 0.0070 0.0286 0.0425 0.0368

?

[University Glasgow]

This is what I want thank you.

[University Glasgow]

I have revised the code but I'm still receiving an error: Arrays have incompatible sizes for this operation

Error in Case1_Ijuptilk (line 28) r=q-(1-alpha).*tan(q)+ C.*tan(q);

%% initialization

clear

clc

close all

%% Model parameters

global k1 eta1 alpha3 gamma1 d N xi C q qstar Theta OddAsymptotes m ...

qq r1 r2 alpha aa b Efun F B Ainv r v theta_sol Vfun v_sol p ModifiedOddAsym ...

qstarTemporary theta_init z AllAsymptotes gg

k1 = 6*10^(-12); % elastic constant

eta1 = 0.0240; % viscosity

xi = 0.01; % activity parameter

alpha3 = -0.001104; % viscosity

gamma1 = 0.1093; % viscosity

Theta = 0.0001;

theta_init = Theta.*sin(pi*z/d);

d = 0.0002;

N = 20;

% Simplified parameters

alpha=1-alpha3^2/gamma1/eta1;

C = (alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);

% Plotting r(q)

q=0:0.01:(5*pi);

r=q-(1-alpha).*tan(q)+ C.*tan(q);

figure

plot(q,r)

axis([0 5*pi -20 20])

% Assign g for q and plot r, and compute the special asymptote

for q=1:N

qstar(q) = 1./(alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);

end

% Compute Odd asymptote

for j=1:N

OddAsymptotes(j)= 1/2*(2*j + 1)*pi;

end

ModifiedOddAsym = [abs(OddAsymptotes) abs(qstar)];

qstarTemporary = min(ModifiedOddAsym);

% Compute the rest of the Roots

for j=1:N

OddAsymptotes(j)= 1/2.*(2*j + 1)*pi;

end

% Compute all symptotes

r = @(q, C) q-(1-alpha).*tan(q)+ C.*tan(q);

AllAsymptotes = [sort(OddAsymptotes) sort(qstar)];

q0 =0;

for q = 1:AllAsymptotes(1)

r1(q) = fzero(@(q)r(q,C), q0);

end

for i=1:N

for q = 2:N

r2(q) = fzero(@(q)r(q,C), q0);

end

end

% Compute all roots

gg = [r1 r2];

% Remove the repeated roots if any & redefine n

qq = unique(gg);

m = length(qq);

% Compute the time constants tan_n

for i=1:m

p(i) = (gamma1*alpha)./(4*k1*qq(i)^2/d^2-alpha3*xi/eta1);

end

% Compute theta_n from initial conditions

Efun = @(z, i, d)cos(qq(i)-2*qq(i).*z/d)-cos(qq(i));

% Compute integral coefficients for off-diagonal elements θ_n matrix

for i=1:m

for j= i+1:m

b(i, j) = integral(@(x)Efun(x, i, d).*@(x)Efun(x, j, d), 0, d);

b(j, i) = b(i, j);

aa(i, j)= b(i, j);

aa(j, i)= b(j, i);

end

end

% Calculate integral coefficients for diagonal elements in theta_n matrix

for i=1:m

aa(i, i)= integral((@(x)Efun(x, i, d)).^2, 0, d);

end

% Calculate integrals for RHS vector

for i=1:m

F(i) = integral(theta_init.*@(x)Efun(i, x, d), 0, d);

end

% Create matrix A and RHS vector B

for i=1:m

for j= 1:m

aa(i, j);

end

end

for i=1:m

B = F(i);

end

% Calculate inverse of A and solve A*x=B

Ainv = inv(A);

r = B.*Ainv;

% Define Theta(z,t)

for i=1:m

theta_sol(i) = sum(r(i).*@(x)Efun(i, x, d).*exp(-t/p(i)));

end

% Compute v_n for n times constant

for i=1:m

v(i) = (-2*k1*alpha3*qq(i))*(1/(d^2*eta1*alpha*gamma1))+ alpha3^2*xi/(2*(eta1)^2*qq(i)*alpha*gamma1)+xi/(2*eta1*qq(i));

end

% Compute v(z,t) from initial conditions

for i=1:m

Vfun(i) = d*sin(qq(i)-2*qq(i)*z/d)+(2*z-d)*sin(qq(i));

end

% Define v(z,t)

for i=1:m

v_sol = sum(v(i).*r(i)*Vfun(i).*exp(-t/p(i)));

end

[Torsten]

z is undefined.

q is defined after you refer to it in the definition of C. This will throw another error.

[University Glasgow]

Thank you for your support. This revised version. I want to use sum function to compute for theta_sol on line 113 and vsol on line 125. I expect to get vectors then plot them but theta_sol and vsol is 0.

%% initialization

clear

clc

close all

%% Model parameters

global k1 eta1 alpha3 gamma1 d N xi C q qstar Theta OddAsymptotes m ...

qq alpha aa b Efun F B r p ModifiedOddAsym qstarTemporary ...

AllAsymptotes gg Efun2 Efun1 v theta_sol Vfun Efun3 v_sol

k1 = 6*10^(-12);            %  elastic constant
eta1 = 0.0240;              %  viscosity
xi = 0.01;                  %  activity parameter
alpha3 = -0.001104;         %  viscosity
gamma1 = 0.1093;            %  viscosity
Theta = 0.0001;
d = 0.0002;
N = 20; 
% Simplified parameters
alpha=1-alpha3^2/gamma1/eta1;
C = (alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);

Unrecognized function or variable 'q'.

%% Plotting r(q)
t = [0:0.1:10];
nt=length(t);
q=0:0.01:(5*pi);
r=q-(1-alpha).*tan(q)+ C.*tan(q);
figure
plot(q,r)
axis([0 5*pi -20 20])
% Assign g for q and plot r, and compute the special asymptote
r=q-(1-alpha).*tan(q)+ C.*tan(q);
for q=1:N
    qstar(q) = 1./(alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);
end
% Compute Odd asymptote
for j=1:N
    OddAsymptotes(j)= 1/2*(2*j + 1)*pi;
end
ModifiedOddAsym = [abs(OddAsymptotes) abs(qstar)];
qstarTemporary = min(ModifiedOddAsym);
% Compute the rest of the Roots
for j=1:N
    OddAsymptotes(j)= 1/2.*(2*j + 1)*pi;
end
% Compute all symptotes
AllAsymptotes = [sort(OddAsymptotes) sort(qstar)];
gg = AllAsymptotes;
% Remove the repeated roots if any & redefine n 
qq = unique(gg);
m = length(qq);
% Compute the time constants tan_n
for i=1:m
    p(i) = (gamma1*alpha)./(4*k1*qq(i)^2/d^2-alpha3*xi/eta1);
end 
% Compute theta_n from initial conditions 
Efun = @(z, i, d)cos(qq(i)-2*qq(i).*z/d)-cos(qq(i));
Efun2 = @(z, j, d)cos(qq(i)-2*qq(j).*z/d)-cos(qq(j));

% Compute integral coefficients for off-diagonal elements θ_n matrix
b = zeros(m, m); 
for i=1:m
    for j= i+1:m
        b(i, j) = integral(@(z)(Efun(z, i, d).*Efun2(z, j, d)), 0, d);
        b(j, i) = b(i, j); 
        aa(i, j) = b(i, j); 
        aa(j, i) = b(j, i);
    end
end
% Calculate integral coefficients for diagonal elements in theta_n matrix
for i=1:m
    aa(i, i)= integral((@(z)Efun(z, i, d).^2), 0, d);
end
% Calculate integrals for RHS vector
Efun1 = @(z, i, d, Theta) Theta.*sin(pi*z/d).*cos(qq(i)-2*qq(i).*z/d)-cos(qq(i));
for i=1:m 
    F(i) = integral(@(z)Efun1(z, i, d, Theta), 0, d);
end  
% Create matrix A and RHS vector BN_mat = N-2;
for i=1:m
    for j= 1:m
        A(i, j) = aa(i, j);
        B(i) = F(i);
    end
end
% Calculate inverse of A and solve A*x=B
Ainv = inv(A);
r = B.*Ainv;
% Compute theta(z,t) from initial conditionsDefine Theta(z,t)
for i =1:m 
    for z=0:d
        Efun3(i) = cos(qq(i)-2*qq(i).*z/d)-cos(qq(i));
    end
end

for i=1:m
    theta_sol(i) = sum(r(i).*Efun3(i).*exp(-t/p(i)));
end
% Compute v_n for n times constant
for i=1:m 
    v(i) = (-2*k1*alpha3*qq(i))*(1/(d^2*eta1*alpha*gamma1))+ alpha3^2*xi/(2*(eta1)^2*qq(i)*alpha*gamma1)+xi/(2*eta1*qq(i));
end
% Compute v(z,t) from initial conditions and define v(z,t)
for i=1:m
    for z=0:d
        Vfun(i) = d*sin(qq(i)-2*qq(i)*z/d)+(2*z-d)*sin(qq(i));
        v_sol = sum(v(i).*r(i)*Vfun(i).*exp(-t/p(i)));
    end 
end 

[Torsten]

for z=0:d

Since d = 0.0002, this loop is empty except 0. Thus Efun3 and v_sol will be 0.

[University Glasgow]

But i want to solve for the theta_v and v_sol for various times for z=0:d. I can help out please.

[Walter Roberson]

Remember that the defined increment is 1, so z=0:d is:

• empty if d is empty
• empty if d(1)<0
• [0] if d(1) < 1
• [0 1] if 1 <= d(1) < 2
• [0 1 2] if 2 <= d(1) < 3

and so on, integers only.

WIth your d = 0.0002, the d(1)<1 case applies, and you would get only [0]

[University Glasgow]

This is the code I tried to translate to matlab. This is actually maple code but maple is limited at this point. I tried to use codegenerator in maple to translate the code but it didn't work. I'm entirely new to matlab. I started using matlab two weeks ago. Please can someone help me to fix the code using my code above.

doCalc:= proc( xi )

# Import Packages
uses ArrayTools, Student:-Calculus1, LinearAlgebra,
ListTools, RootFinding, plots, ListTools:
local gamma__1:= .1093,
alpha__3:= -0.1104e-2,
k__1:= 6*10^(-12),
d:= 0.2e-3,
theta0:= 0.001,
eta__1:= 0.240e-1,
alpha:= 1-alpha__3^2/(gamma__1*eta__1), 
c:= alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha__3*xi/eta__1)),
theta_init:= theta0*sin(Pi*z/d),
n:= 30,
g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,
qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,
soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,
AllAsymptotes, p, Efun, b, aa, F, A, B, Ainv, r, theta_sol, v, Vfun, v_sol,minp,nstar;

# Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes

g:= q-(1-alpha)*tan(q)+ c*tan(q):
f:= subs(q = x+I*y, g):
b1:= evalc(Re(f)) = 0: 
b2:= evalc(Im(f)) = 0:
qstar:= (fsolve(1/c = 0, q = 0 .. infinity)):
OddAsymptotes:= Vector[row]([seq(evalf(1/2*(2*j + 1)*Pi), j = 0 .. n)]);

# Compute Odd asymptote

ModifiedOddAsym:= abs(`-`~(OddAsymptotes, qstar));
qstarTemporary:= min(ModifiedOddAsym);
indexOfqstar2:= SearchAll(qstarTemporary, ModifiedOddAsym);
qstar2:= OddAsymptotes(indexOfqstar2);

# Compute complex roots

AreThereComplexRoots:= type(true, 'truefalse');
try
    #print({min(qstar2, qstar), max(qstar2, qstar)});
    #print({eval(subs({x=min(qstar2, qstar)+0.1,y=0},b1))});
    soln1:= fsolve({b1, b2}, {x = min(qstar2, qstar)+0.1 .. max(qstar2, qstar)-0.1, y = 0.1 .. 100}); 
    soln2:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = -infinity .. 0}); 
    #print(soln1);
    qcomplex1:= subs(soln1, x+I*y); 
    qcomplex2:= subs(soln2, x+I*y);
catch:
    AreThereComplexRoots:= type(FAIL, 'truefalse');
end try;
    
    # Compute the rest of the Roots
    OddAsymptotes:= Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]); 
    AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar])); 
    if AreThereComplexRoots
        then gg:= [ qcomplex1, qcomplex2,op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),
        seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
        elif not AreThereComplexRoots 
        then gg:= [op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] ..                           AllAsymptotes[i+1])), i = 1 .. n)];
    end if:
        
        # Remove the repeated roots if any & Redefine n
        
        qq:= MakeUnique(gg):
        m:= numelems(qq):
        
        ## Compute the 
        `τ_n`;
        time constants
        
        for i to m do 
            p[i] := gamma__1*alpha/(4*k__1*qq[i]^2/d^2-alpha__3*xi/eta__1):
        end do:
        
        ## Compute θ_n from initial conditions
        
        for i to m do 
            Efun[i] := cos(qq[i]-2*qq[i]*z/d)-cos(qq[i]):
        end do:
        
        ## Compute integral coefficients for off-diagonal elements θ_n matrix
        
        printlevel := 2:
        for i to m do 
            for j from i+1 to m do 
                b[i, j] := int(Efun[i]*Efun[j], z = 0 .. d): 
                b[j, i] := b[i, j]: 
                aa[i, j] := b[i, j]: 
                aa[j, i] := b[j, i]:
            end do :
        end do:
        
        ## Calculate integral coefficients for diagonal elements in theta_n matrix
        
        for i to m do 
            aa[i, i] := int(Efun[i]^2, z = 0 .. d): 
        end do:
        
        ## Calculate integrals for RHS vector
        
        for i to m do 
            F[i] := int(theta_init*Efun[i], z = 0 .. d):
        end do:
        
        ## Create matrix A and RHS vector B
        
        A := Matrix([seq([seq(aa[i, j], i = 1 .. m)], j = 1 .. m)]): 
        B := Vector([seq(F[i], i = 1 .. m)]):
        
        ## Calculate inverse of A and solve A*x=B
        
        Ainv := 1/A:
        r := MatrixVectorMultiply(Ainv, B):
        
        ## Define Theta(z,t)
        theta_sol := add(r[i]*Efun[i]*exp(-t/p[i]), i = 1 .. m):
        
        ## Compute v_n for n times constant
        
        for i to m do 
            v[i] := (-2*k__1*alpha__3*qq[i])*(1/(d^2*eta__1*alpha*gamma__1))+ alpha__3^2*xi/(2*(eta__1)^2*qq[i]*alpha*gamma__1)+xi/(2*eta__1*qq[i]):
        end do:
        
        ## Compute v(z,t) from initial conditions
        for i to m do 
            Vfun[i] := d*sin(qq[i]-2*qq[i]*z/d)+(2*z-d)*sin(qq[i]):
        end do:
        
        ## Define v(z,t)
        v_sol := add(v[i]*r[i]*Vfun[i]*exp(-t/p[i]), i = 1 .. m):
        
        ## 
        minp:=min( seq( Re(p[i]), i=1..m) ):
        member(min(seq( Re(p[i]), i=1..m)),[seq( Re(p[i]), i=1..m)],'nstar'):
        
        ## Return all the plots
        return theta_init, theta_sol, v_sol, minp, eval(p), nstar, p[nstar], g, qq, qstar, AreThereComplexRoots;
    end proc:

[Steven Lord]

You've posted a lot of code. But let's take a step back for a second. Can you explain in words not code the problem you're trying to solve? What do you know and what are you trying to find or compute? Knowing the context may help us offer better guidance in how to adjust the code.

[University Glasgow]

First of all, I have a problem with computing the root of the function r. It seems that maple do skip some roots of the function r that's why I have change to matlab. I want to compute the roots of r for q=(2*i-2)*Pi/2..(2*i+1)*Pi/2) and assign the results to gg, if there is any repeated, i used unique() to remove the repeated roots and finally assign the result to qq and use qq to compute p(i), theta_sol and v_sol. I don't know how to find roots with matlab. I tried vpasolve() but i need syms which didn't work in my case.

Finally, I want to plot theta_sol and v_sol for various times at z=0:d.

%% initialization
clear 
clc
close all

%% Model parameters
global k1 eta1 alpha3 gamma1 N xi C q qstar Theta OddAsymptotes m ...
    qq alpha aa b Efun F B  r p ModifiedOddAsym qstarTemporary ...
    AllAsymptotes gg Efun2 Efun1 v theta_sol Vfun Efun3 v_sol d z
k1 = 6*10^(-12);            %  elastic constant
eta1 = 0.0240;              %  viscosity
xi = 0.01;                  %  activity parameter
alpha3 = -0.001104;         %  viscosity
gamma1 = 0.1093;            %  viscosity
Theta = 0.0001;
d = 0.0002;
N = 20; 
% Simplified parameters
alpha=1-alpha3^2/gamma1/eta1;
C = (alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1);
%% Plotting r(q)
t = [0:0.1:10];
nt=length(t);
q=0:0.01:(5*pi);
r=q-(1-alpha).*tan(q)+ C.*tan(q);
figure
plot(q,r)
axis([0 5*pi -20 20])

% Remove the repeated roots if any & redefine n 
qq = unique(gg);
m = length(qq);

% Compute the time constants tan_n
for i=1:m
    p(i) = (gamma1*alpha)./(4*k1*qq(i)^2/d^2-alpha3*xi/eta1);
end 
% Compute theta_n from initial conditions 
Efun = @(z, i, d)cos(qq(i)-2*qq(i).*z/d)-cos(qq(i));
Efun2 = @(z, j, d)cos(qq(i)-2*qq(j).*z/d)-cos(qq(j));

% Compute integral coefficients for off-diagonal elements θ_n matrix
b = zeros(m, m); 
for i=1:m
    for j= i+1:m
        b(i, j) = integral(@(z)(Efun(z, i, d).*Efun2(z, j, d)), 0, d);
        b(j, i) = b(i, j); 
        aa(i, j) = b(i, j); 
        aa(j, i) = b(j, i);
    end
end
% Calculate integral coefficients for diagonal elements in theta_n matrix
for i=1:m
    aa(i, i)= integral((@(z)Efun(z, i, d).^2), 0, d);
end
% Calculate integrals for RHS vector
Efun1 = @(z, i, d, Theta) Theta.*sin(pi*z/d).*cos(qq(i)-2*qq(i).*z/d)-cos(qq(i));
for i=1:m 
    F(i) = integral(@(z)Efun1(z, i, d, Theta), 0, d);
end  
% Create matrix A and RHS vector BN_mat = N-2;
for i=1:m
    for j= 1:m
        A(i, j) = aa(i, j);
        B(i) = F(i);
    end
end
% Calculate inverse of A and solve A*x=B
Ainv = inv(A);
b = Ainv\B';
% Compute theta(z,t) from initial conditionsDefine Theta(z,t)              
for i=1:m
    z = 0:d;
    Efun3(i) = cos(qq(i)-2*qq(i).*z/d)-cos(qq(i));
    theta_sol(i) = sum(b(i).*Efun3(i).*exp(-t/p(i)));
end
% Compute v_n for n times constant
for i=1:m 
    v(i) = (-2*k1*alpha3*qq(i))*(1/(d^2*eta1*alpha*gamma1))+ alpha3^2*xi/(2*(eta1)^2*qq(i)*alpha*gamma1)+xi/(2*eta1*qq(i));
end
% Compute v(z,t) from initial conditions and define v(z,t)
for i=1:m
    Vfun(i) = d*sin(qq(i)-2*qq(i)*z/d)+(2*z-d)*sin(qq(i));
    v_sol = sum(v(i).*r(i)*Vfun(i).*exp(-t/p(i)));
end