I am solving nonlinear PDEs and receive this error:Error using daeic12 (line 72) This DAE appears to be of index greater than 1.

The lcode_nondim.m defines the right-hand side of the dimensionless PDEs for theta (director orientation) and v (flow velocity). It takes the current state and parameters, computes the nonlinear coefficients and finite difference.

Max_Non_linear_Nondim_func.m sets up the problem and solves the dimensionless PDEs defined in lcode_nondim.m. It initializes parameters, constructs the mass matrix needed for the DAE solver ode23t, sets initial conditions for theta and v, and defines the nondimensional grid in space and time. Then it calls lcode_nondim through ode23t to evolve the system in time. After solving, it reconstructs the full solutions including boundary values and extracts the values of theta at specific positions (middle, quarter, and three-quarter points) for later analysis or plotting.

Results_rho0.m is the main script that initializes the activity parameter xi and the magnetic field parameter rho. It calls Max_Non_linear_Nondim_func to run the simulation for these parameters. After the simulation, it plots theta versus xi at different positions in the layer (middle, quarter, and three-quarter) to visualize how the director orientation evolves with the activity. The script sets up figure properties such as grid, labels, legend, and line styles for clarity.

The only file to run is the Result_rho0.m

%% initialization
clear 
clc
close all

%% Activity parameter values
xi = 0;  % dimensionless activity

%% Magnetic field parameter (dimensionless rho)
rho = 0.0;  % corresponds to F/F_c

%% Run dimensionless simulation
[THETA_dimless, V_dimless, theta_middle_tmax_dimless, theta_middle_tmax41_dimless, theta_middle_tmax_3d4_dimless, t, z, h, N, d] = ...
    Max_Non_linear_Nondim_func(xi, rho);

Error using daeic12 (line 72)
This DAE appears to be of index greater than 1.

Error in ode23t (line 278)
[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(ode,odeArgs,t,ICtype,Mt,y,yp0,f0,...
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in solution>Max_Non_linear_Nondim_func (line 84)
[t, y] = ode23t(@(t,y) lcode_nondim(t, y, xi, rho, N, Phi_b, h, kappa, eta2hat, eta12hat, gamma2hat, alpha), tspan, u0, options);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

%% Plot theta at different positions vs xi
figure; 
hold on; box on; grid on;

plot(xi, theta_middle_tmax_dimless, 'r-', 'LineWidth', 2);
plot(xi, theta_middle_tmax41_dimless, 'b--', 'LineWidth', 2);
plot(xi, theta_middle_tmax_3d4_dimless, 'g-.', 'LineWidth', 2);

xlabel('\xi (dimensionless activity)');
ylabel('\theta (dimensionless)');
title('\theta at different layer positions vs \xi');
legend('Middle (z=1/2)','Quarter (z=1/4)','3/4 (z=3/4)','Location','best');

set(gca,'FontSize',14);

function [THETA, V, theta_middle_tmax, theta_middle_tmax41, ...
    theta_middle_tmax_3d4, t, z, h, N, d] = Max_Non_linear_Nondim_func(xi, rho)

%% ---------------- Nondimensional parameters ----------------
gamma1 = 0.1093;
gamma2 = -0.1121;
eta1   = 0.0240;
eta2   = 0.1361;
eta12  = -0.0181;
K1     = 6e-12;
K3     = 7.5e-12;

alpha      = (gamma1 + gamma2)^2 / (4*gamma1*eta1);
eta2hat    = eta2 / eta1;
eta12hat   = eta12 / eta1;
gamma2hat  = gamma2 / (gamma1 + gamma2);
kappa      = K3 / K1;

d = 0.0002;   % Layer thickness
N = 200;      % interior points
h = 1/N;      % nondimensional spacing

z = linspace(0,1,N+1);  % nondimensional z
tspan = 0:2000;         % nondimensional time

%% ---------------- Initial condition ----------------
Theta0 = 0.0001;
theta0 = Theta0*(sin(pi*z)+0.5*sin(2*pi*z));
v0     = zeros(1,N+1);

theta0_int = theta0(2:N);
v0_int     = v0(2:N);
u0 = [theta0_int'; v0_int'];

%% ---------------- Mass matrix for ode23t ----------------
M1 = eye(N-1);
M2 = zeros(N-1);
M = [M1 M2; M2 M2];

%% ---------------- Boundary conditions ----------------
Phi_b = 0; % theta at boundaries
v_b = 0;   % v at boundaries

%% ---------------- Solve the PDEs ----------------
counter = 0;
Xis = xi;  % can handle vector of xi values
for xiidx = 1:length(Xis)
    xi = Xis(xiidx);
    
    options = odeset('Mass', M, 'RelTol', 1e-4, 'AbsTol', 1e-6);
    
    % Solve using ode23t with dimensionless PDEs
    [t, y] = ode23t(@(t,y) lcode_nondim(t, y, xi, rho, N, Phi_b, h, kappa, eta2hat, eta12hat, gamma2hat, alpha), tspan, u0, options);
    
    % Reconstruct full theta and v (including boundaries)
    theta = [Phi_b*ones(length(t),1) y(:,1:N-1) Phi_b*ones(length(t),1)];
    v     = [v_b*ones(length(t),1) y(:,N:2*N-2) v_b*ones(length(t),1)];

    % Store solutions
    counter = counter + 1;
    THETA(counter,:,:) = theta;
    V(counter,:,:) = v;

    % theta at z = d/2
    theta_middle(counter,:,:) = theta(:, N/2+1);
    theta_middle_tmax(counter) = theta_middle(counter,end);

    % theta at z = d/4
    theta_middle41(counter,:,:) = theta(:, N/4+1);
    theta_middle_tmax41(counter) = theta_middle41(counter,end);

    % theta at z = 3d/4
    theta_middle_3d(counter,:,:) = theta(:, 3*N/4+1);
    theta_middle_tmax_3d4(counter) = theta_middle_3d(counter,end);
end

end

%% Functions
% Right hand side of the ODEs for dimensionless equations
function rhsode = lcode_nondim(t, y, xi, rho, N, Phi_b, h, kappa, eta2_hat, eta12_hat, gamma2_hat, alpha)

% initialize theta and v
theta = y(1:(N-1));
v     = y(N:(2*N-2));

% --- Simplified parameters ---
gbar     = cos(theta).^2 + eta2_hat*sin(theta).^2 + eta12_hat*sin(theta).^2.*cos(theta).^2;
gbar_der = -2*cos(theta).*sin(theta) + 2*eta2_hat*cos(theta).*sin(theta) + 2*eta12_hat*cos(theta).*sin(theta).*(cos(theta).^2 - sin(theta).^2);

mbar     = 1 - 2*gamma2_hat*sin(theta).^2;
mbar_der = -4*gamma2_hat*sin(theta).*cos(theta);

f1       = cos(theta).^2 + kappa*sin(theta).^2;
f2       = (kappa-1)*sin(theta).*cos(theta);

f1_der   = -2*cos(theta).*sin(theta) + 2*kappa*sin(theta).*cos(theta);
f2_der   = (kappa-1)*(cos(theta).^2 - sin(theta).^2);

% --- Derivatives for v-equation nonlinear terms ---
Kappa1 = gbar_der;  % derivative of gbar
% Kappa2 = mbar_der;  % derivative of mbar
% Kappa3 = f1_der;    % derivative of f1 in theta_tau
% Kappa4 = f2_der;    % derivative of f2 in theta_tau

% Coefficients
A1 = f1;
B1 = f2;
C1 = mbar;

A2 = alpha.*mbar_der.*A1;
B2 = alpha.*mbar_der.*B1;
C2 = alpha.*mbar_der.*C1;
D2 = alpha.*mbar_der.*(pi^2*rho^2);
A3 = alpha.*f1_der.*C1;
B3 = A1.*C1;
B4 = B1.*C1;
C3 = f2_der.*C1;
C4 = C1.*C1;
E1 = (pi^2*rho^2).*C1.*cos(2*theta);
E2 = 4*pi^2*xi*cos(2*theta);

% Initialize RHS
%rhsode = zeros(2*(N-1),1);

%% --- Theta equations ---
% Left boundary
rhsode(1,1) = (A1(1)/(h^2)).*(theta(2)-2*theta(1)+Phi_b) ...
             + B1(1).*((theta(2)-Phi_b)./(2*h))^2 ...
             - (C1(1)/(2*h)).*(v(2)-0) ...
             + pi^2*rho^2.*sin(theta(1))*cos(theta(1));

% Interior points
for i = 2:N-2
    rhsode(i,1) = (A1(i)./h^2)*(theta(i+1)-2*theta(i)+theta(i-1)) ...
                 + B1(i).*((theta(i+1)-theta(i-1))./(2*h))^2 ...
                 - (C1(i)/(2*h)).*(v(i+1)-v(i-1)) ...
                 + pi^2*rho^2.*sin(theta(i))*cos(theta(i));
end

% Right boundary
rhsode(N-1,1) = (A1(N-1)/h^2)*(Phi_b-2*theta(N-1)+theta(N-2)) ...
               + B1(N-1).*((Phi_b-theta(N-2))./(2*h))^2 ...
               - (C1(N-1)./(2*h)).*(0-v(N-2)) ...
               + pi^2*rho^2.*sin(theta(N-1))*cos(theta(N-1));
%% --- v equations ---
% Left boundary points
rhsode(N,1) = (B3(1)/(h^3)).*(-Phi_b+3*theta(1)-3*theta(2)+theta(3)) ...
             + ((gbar(1)-C4(1))./h^2).*(v(2)-2*v(1)+0) ...
             + (A2(1)+A3(1)+2*B4(1))*((theta(2)-Phi_b)/(2*h))*((theta(2)-2*theta(1)+Phi_b)/h^2) ...
             + ((Kappa1(1)-2*C2(1))*(v(2)-0)/(2*h) + D2(1).*sin(theta(1))*cos(theta(1)) + E1(1)+E2(1))*((theta(2)-Phi_b)/(2*h)) ...
             + (B2(1)+C3(1))*((theta(2)-Phi_b)/(2*h))^3;

rhsode(N+1,1) = (B3(2)/(2*h.^3)).*(-Phi_b+2*theta(1)-2*theta(3)+theta(4)) ...
               + ((gbar(2)-C4(2))./h^2).*(v(3)-2*v(2)+v(1)) ...
               + (A2(2)+A3(2)+2*B4(2))*((theta(3)-theta(1))/(2*h))*((theta(3)-2*theta(2)+theta(1))/h^2) ...
               + ((Kappa1(2)-2*C2(2))*(v(3)-v(1))/(2*h) + D2(2).*sin(theta(2))*cos(theta(2)) + E1(2)+E2(2))*((theta(3)-theta(1))/(2*h)) ...
               + (B2(2)+C3(2))*((theta(3)-theta(1))/(2*h))^3;



% v equations [REMEMBER theta the RHS index for the v equations is N-1 more than the indices of the theta and v variables in this function]
% for the two positons to the rigbar ht of the left hand boundary z=0
rhsode(N,1)     =  (B3(1)/(h.^3)).*(-Phi_b + 3*theta(1) - 3*theta(2) + theta(3)) ...
                   + ((gbar (1)-C4(1))./(h^2)).*(v(2) -2*v(1)+ 0) ...
                   + ( A2(1)  +   A3(1)   +2*B4(1) )*((theta(2)-Phi_b)/(2*h))*(theta(2)-2*theta(1)+ Phi_b)/(h^2)...
                   + ( (Kappa1(1)-2*C2(1))*(1/(2*h))*(v(2)-0)+ D2(1).*sin(theta(1))*cos(theta(1))   + E1(1) + E2(1))* (1/(2*h))*(theta(2)-Phi_b)...
                   + ( B2(1)+ C3(1))*(( 1/(2*h)).*(theta(2)-Phi_b)).^3;

rhsode(N+1,1)   = (B3(2)/(2*h.^3)).*(-Phi_b + 2*theta(1) - 2*theta(3) + theta(4)) ...
                   + ( (gbar (2)-C4(2))/(h^2)).*(v(3) -2*v(2)+ v(1) ) ...
                   + ( A2(2)  +A3(2)  +2*B4(2) )*((theta(3)-theta(1))/(2*h)).*(theta(3)-2*theta(2)+ theta(1))/(h^2)...
                  + ( (Kappa1(2) -    2*C2(2))*(1/(2*h))*(v(3)-v(1))   + D2(2).*sin(theta(2))*cos(theta(2))   + E1(2) + E2(2) )*(1/(2*h))*(theta(3)-theta(1))...
                  + ( B2(2)   + C3(2))*( (1/(2*h))*(theta(3)-theta(1))).^3;

% for all other internal positions 0<z<d
for i=(N+2):(2*N-4)
    rhsode(i,1) = (B3(i-(N-1))/(2*h.^3))*(-theta(i-2-(N-1)) + 2*theta(i-1-(N-1)) - 2*theta(i+1-(N-1))+ theta(i+2-(N-1)))         + ( (gbar (i-(N-1))-C4(i-(N-1)))/(h.^2))*( v(i+1-(N-1)) - 2*v(i-(N-1)) + v(i-1-(N-1)) )                             + ( A2(i-(N-1))   +A3(i-(N-1)) +2*B4(i-(N-1)))*((theta(i+1-(N-1))-theta(i-1-(N-1)))/(2*h)).*(theta(i+1-(N-1)) - 2*theta(i-(N-1)) + theta(i-1-(N-1)))/(h^2)...
                  + ( (Kappa1(i-(N-1))  -   2*C2(i-(N-1)))*(1/(2*h))*( v(i+1-(N-1))-v(i-1-(N-1)))   + D2(i-(N-1))*sin(theta(i-(N-1)))*cos(theta(i-(N-1)))   + E1(i-(N-1)) + E2(i-(N-1)))*(1/(2*h))*( theta(i+1-(N-1))-theta(i-1-(N-1)))         + (B2(i-(N-1))   + C3(i-(N-1)))*((1/(2*h))*( theta(i+1-(N-1))-theta(i-1-(N-1)))).^3;

end
% for the two positons to the left of the rigbar ht hand boundary z=d
rhsode(2*N-3,1) = (B3(N-2)/(2*h^3))*( -theta(N-4) + 2*theta(N-3) - 2*theta(N-1) + Phi_b )                                + ((gbar (N-2)- C4(N-2))/h.^2)*(v(N-1) -2*v(N-2) + v(N-3))                  + ( A2(N-2)  +A3(N-2)  +2*B4(N-2) )*((theta(N-1)-theta(N-3))/(2*h)).*(theta(N-1) -2*theta(N-2) + theta(N-3))/(h^2)...
                  + ((Kappa1(N-2) -   2*C2(N-2))*(1/(2*h))*( v(N-1)-v(N-3) ) + D2(N-2)*sin(theta(N-2)).*cos(theta(N-2))  + E1(N-2)  + E2(N-2) )*(1/(2*h))*(theta(N-1)-theta(N-3))                + ( B2(N-2)  + C3(N-2) )*((1/(2*h))*(theta(N-1)-theta(N-3))).^3;

rhsode(2*N-2,1) = (B3(N-1)/(h^3)).*( -theta(N-3) + 3*theta(N-2) - 3*theta(N-1) + Phi_b )                                 + (( gbar (N-1)-C4(N-1))./(h.^2)).*( 0 -2*v(N-1) + v(N-2))                  + (A2(N-1) +A3(N-1) +2*B4(N-1) )* (( Phi_b-theta(N-2))/(2*h)).*( Phi_b -2*theta(N-1) + theta(N-2))/(h^2)...
                  + ( (Kappa1(N-1) - 2*C2(N-1))*(1/(2*h))*(0-v(N-2))   + D2(N-1).*sin(theta(N-1)).*cos(theta(N-1))       + E1(N-1) + E2(N-1))*(1/(2*h))*(  Phi_b-theta(N-2))                     + ( B2(N-1)  + C3(N-1) )*((1/(2*h))*( Phi_b-theta(N-2))).^3;
                            
end 

[dpb]

AI-generated:

The error "This DAE appears to be of index greater than 1" with ode23 is caused by the DAE system not being index-1, which is a requirement for the solver, often due to inconsistent equations, improper scaling, or an incorrect initial guess. To fix this, you must ensure your initial conditions (y0) are consistent with the DAEs, check for contradictions in your equations, and improve numerical stability by scaling your variables.

Common causes and solutions

• Contradictory or inconsistent equations: Your algebraic equations may contain contradictions that make them unsolvable.
• Solution: Carefully review your equations for any algebraic contradictions. For example, two equations for the same variable might lead to an impossible relationship.
• Inconsistent initial conditions: Your initial conditions (y0) might not satisfy the algebraic equations of your DAE system.
• Solution: Ensure your y0 values are consistent with the algebraic constraints of your system. You may need to solve for a consistent initial condition before starting the simulation.
• Poor scaling: Variables with vastly different magnitudes can lead to numerical issues where small terms are neglected, effectively increasing the DAE index.
• Solution: Rescale your variables so they have a similar magnitude. This can involve dividing by a constant or using odeset to adjust absolute and relative tolerances ('RelTol', 'AbsTol').
• Numerical singularity: The system's Jacobian can become singular for certain parameters, leading to numerical instability that the solver interprets as an index>1 DAE.
• Solution: Try changing the parameters. A change in a parameter's value from 10E-10 to 10E-12, for example, could trigger the error.
• Incorrect use of solvers and options: Some solvers and options interact in unexpected ways. For instance, setting the mass matrix M to 'sparse' can change how the DAE is treated and may cause issues.
• Solution: Experiment with different solvers, like ode15s, and ensure your options are set correctly. For example, using ode15s with 'MassSingular', 'yes' might help if the issue is related to the mass matrix.

Without the actual system trying to be solved written out, tough for anybody else to have much chance.

[University]

Thank you for your invaluable contribution. I will check your suggestions.

[Torsten]

@University

It would help if you included the PDEs you are trying to solve (equations, boundary conditions, initial conditions) and the discretization you use in your code.

[University]

Here is it

[Torsten]

Often you get this error message because the integrator cannot determine consistent initial values for the algebraic equations.

You define the first (N-1) equations as ODEs for "theta" and the last (N-1) equations as algebraic equations. Since "theta" and "v" are the unknowns, the last (N-1) algebraic equations must be uniquely solvable for v(1),...,v(N-1). Before calling ode23t, you should test this. So fix "theta" to the initial values you supply (theta0 = Theta0*(sin(pi*z)+0.5*sin(2*pi*z))) and use a solver for nonlinear algebraic equations (e.g. "fsolve") to determine v(1),...,v(N-1) from the last (N-1) algebraic equations. If "fsolve" doesn't succeed, "ode23t" won't succeed either. If "fsolve" succeeds, supply the values for v(1),...,v(N-1) as initial values for v for "ode23t" instead of v == 0 as you do at the moment.

[University]

Okay, thank you. I will let you know what I get out of it.

[Torsten]

I computed the solution for v0 corresponding to your initial theta0 and get values in the order of 1e11.

So I guess there must be something wrong with your algebraic equations.

%% Activity parameter values
xi = 0;  % dimensionless activity

%% Magnetic field parameter (dimensionless rho)
rho = 0.0;  % corresponds to F/F_c

%% Run dimensionless simulation
Max_Non_linear_Nondim_func(xi, rho);

ans = 199

ans = 199

ans = 199×1

1.0e+11 * 0.0221 0.2009 1.8293 1.8474 1.8655 1.8836 1.9018 1.9199 1.9381 1.9563 1.9745 1.9927 0.2200 2.0110 2.0292

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ans = 199×1

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function Max_Non_linear_Nondim_func(xi, rho)

%% ---------------- Nondimensional parameters ----------------
gamma1 = 0.1093;
gamma2 = -0.1121;
eta1   = 0.0240;
eta2   = 0.1361;
eta12  = -0.0181;
K1     = 6e-12;
K3     = 7.5e-12;

alpha      = (gamma1 + gamma2)^2 / (4*gamma1*eta1);
eta2hat    = eta2 / eta1;
eta12hat   = eta12 / eta1;
gamma2hat  = gamma2 / (gamma1 + gamma2);
kappa      = K3 / K1;

d = 0.0002;   % Layer thickness
N = 200;      % interior points
h = 1/N;      % nondimensional spacing

z = linspace(0,1,N+1);  % nondimensional z

%% ---------------- Initial condition ----------------
Theta0 = 0.0001;
theta0 = Theta0*(sin(pi*z)+0.5*sin(2*pi*z));

theta0_int = theta0(2:N);

%% ---------------- Boundary conditions ----------------
Phi_b = 0; % theta at boundaries

lcode_nondim(0, theta0_int', xi, rho, N, Phi_b, h, kappa, eta2hat, eta12hat, gamma2hat, alpha);

end

%% Functions
% Right hand side of the ODEs for dimensionless equations
function lcode_nondim(t, theta, xi, rho, N, Phi_b, h, kappa, eta2_hat, eta12_hat, gamma2_hat, alpha)

% initialize theta and v
%theta = y(1:(N-1));
%v     = y(N:(2*N-2));
v = sym('v',[N-1,1]);

% --- Simplified parameters ---
gbar     = cos(theta).^2 + eta2_hat*sin(theta).^2 + eta12_hat*sin(theta).^2.*cos(theta).^2;
gbar_der = -2*cos(theta).*sin(theta) + 2*eta2_hat*cos(theta).*sin(theta) + 2*eta12_hat*cos(theta).*sin(theta).*(cos(theta).^2 - sin(theta).^2);

mbar     = 1 - 2*gamma2_hat*sin(theta).^2;
mbar_der = -4*gamma2_hat*sin(theta).*cos(theta);

f1       = cos(theta).^2 + kappa*sin(theta).^2;
f2       = (kappa-1)*sin(theta).*cos(theta);

f1_der   = -2*cos(theta).*sin(theta) + 2*kappa*sin(theta).*cos(theta);
f2_der   = (kappa-1)*(cos(theta).^2 - sin(theta).^2);

% --- Derivatives for v-equation nonlinear terms ---
Kappa1 = gbar_der;  % derivative of gbar
% Kappa2 = mbar_der;  % derivative of mbar
% Kappa3 = f1_der;    % derivative of f1 in theta_tau
% Kappa4 = f2_der;    % derivative of f2 in theta_tau

% Coefficients
A1 = f1;
B1 = f2;
C1 = mbar;

A2 = alpha.*mbar_der.*A1;
B2 = alpha.*mbar_der.*B1;
C2 = alpha.*mbar_der.*C1;
D2 = alpha.*mbar_der.*(pi^2*rho^2);
A3 = alpha.*f1_der.*C1;
B3 = A1.*C1;
B4 = B1.*C1;
C3 = f2_der.*C1;
C4 = C1.*C1;
E1 = (pi^2*rho^2).*C1.*cos(2*theta);
E2 = 4*pi^2*xi*cos(2*theta);

% Initialize RHS
%rhsode = zeros(2*(N-1),1);

%% --- Theta equations ---
% Left boundary
rhsode(1,1) = (A1(1)/(h^2)).*(theta(2)-2*theta(1)+Phi_b) ...
             + B1(1).*((theta(2)-Phi_b)./(2*h))^2 ...
             - (C1(1)/(2*h)).*(v(2)-0) ...
             + pi^2*rho^2.*sin(theta(1))*cos(theta(1));

% Interior points
for i = 2:N-2
    rhsode(i,1) = (A1(i)./h^2)*(theta(i+1)-2*theta(i)+theta(i-1)) ...
                 + B1(i).*((theta(i+1)-theta(i-1))./(2*h))^2 ...
                 - (C1(i)/(2*h)).*(v(i+1)-v(i-1)) ...
                 + pi^2*rho^2.*sin(theta(i))*cos(theta(i));
end

% Right boundary
rhsode(N-1,1) = (A1(N-1)/h^2)*(Phi_b-2*theta(N-1)+theta(N-2)) ...
               + B1(N-1).*((Phi_b-theta(N-2))./(2*h))^2 ...
               - (C1(N-1)./(2*h)).*(0-v(N-2)) ...
               + pi^2*rho^2.*sin(theta(N-1))*cos(theta(N-1));
%% --- v equations ---
% Left boundary points
rhsode(N,1) = (B3(1)/(h^3)).*(-Phi_b+3*theta(1)-3*theta(2)+theta(3)) ...
             + ((gbar(1)-C4(1))./h^2).*(v(2)-2*v(1)+0) ...
             + (A2(1)+A3(1)+2*B4(1))*((theta(2)-Phi_b)/(2*h))*((theta(2)-2*theta(1)+Phi_b)/h^2) ...
             + ((Kappa1(1)-2*C2(1))*(v(2)-0)/(2*h) + D2(1).*sin(theta(1))*cos(theta(1)) + E1(1)+E2(1))*((theta(2)-Phi_b)/(2*h)) ...
             + (B2(1)+C3(1))*((theta(2)-Phi_b)/(2*h))^3;

rhsode(N+1,1) = (B3(2)/(2*h.^3)).*(-Phi_b+2*theta(1)-2*theta(3)+theta(4)) ...
               + ((gbar(2)-C4(2))./h^2).*(v(3)-2*v(2)+v(1)) ...
               + (A2(2)+A3(2)+2*B4(2))*((theta(3)-theta(1))/(2*h))*((theta(3)-2*theta(2)+theta(1))/h^2) ...
               + ((Kappa1(2)-2*C2(2))*(v(3)-v(1))/(2*h) + D2(2).*sin(theta(2))*cos(theta(2)) + E1(2)+E2(2))*((theta(3)-theta(1))/(2*h)) ...
               + (B2(2)+C3(2))*((theta(3)-theta(1))/(2*h))^3;



% v equations [REMEMBER theta the RHS index for the v equations is N-1 more than the indices of the theta and v variables in this function]
% for the two positons to the rigbar ht of the left hand boundary z=0
rhsode(N,1)     =  (B3(1)/(h.^3)).*(-Phi_b + 3*theta(1) - 3*theta(2) + theta(3)) ...
                   + ((gbar (1)-C4(1))./(h^2)).*(v(2) -2*v(1)+ 0) ...
                   + ( A2(1)  +   A3(1)   +2*B4(1) )*((theta(2)-Phi_b)/(2*h))*(theta(2)-2*theta(1)+ Phi_b)/(h^2)...
                   + ( (Kappa1(1)-2*C2(1))*(1/(2*h))*(v(2)-0)+ D2(1).*sin(theta(1))*cos(theta(1))   + E1(1) + E2(1))* (1/(2*h))*(theta(2)-Phi_b)...
                   + ( B2(1)+ C3(1))*(( 1/(2*h)).*(theta(2)-Phi_b)).^3;

rhsode(N+1,1)   = (B3(2)/(2*h.^3)).*(-Phi_b + 2*theta(1) - 2*theta(3) + theta(4)) ...
                   + ( (gbar (2)-C4(2))/(h^2)).*(v(3) -2*v(2)+ v(1) ) ...
                   + ( A2(2)  +A3(2)  +2*B4(2) )*((theta(3)-theta(1))/(2*h)).*(theta(3)-2*theta(2)+ theta(1))/(h^2)...
                  + ( (Kappa1(2) -    2*C2(2))*(1/(2*h))*(v(3)-v(1))   + D2(2).*sin(theta(2))*cos(theta(2))   + E1(2) + E2(2) )*(1/(2*h))*(theta(3)-theta(1))...
                  + ( B2(2)   + C3(2))*( (1/(2*h))*(theta(3)-theta(1))).^3;

% for all other internal positions 0<z<d
for i=(N+2):(2*N-4)
    rhsode(i,1) = (B3(i-(N-1))/(2*h.^3))*(-theta(i-2-(N-1)) + 2*theta(i-1-(N-1)) - 2*theta(i+1-(N-1))+ theta(i+2-(N-1)))         + ( (gbar (i-(N-1))-C4(i-(N-1)))/(h.^2))*( v(i+1-(N-1)) - 2*v(i-(N-1)) + v(i-1-(N-1)) )                             + ( A2(i-(N-1))   +A3(i-(N-1)) +2*B4(i-(N-1)))*((theta(i+1-(N-1))-theta(i-1-(N-1)))/(2*h)).*(theta(i+1-(N-1)) - 2*theta(i-(N-1)) + theta(i-1-(N-1)))/(h^2)...
                  + ( (Kappa1(i-(N-1))  -   2*C2(i-(N-1)))*(1/(2*h))*( v(i+1-(N-1))-v(i-1-(N-1)))   + D2(i-(N-1))*sin(theta(i-(N-1)))*cos(theta(i-(N-1)))   + E1(i-(N-1)) + E2(i-(N-1)))*(1/(2*h))*( theta(i+1-(N-1))-theta(i-1-(N-1)))         + (B2(i-(N-1))   + C3(i-(N-1)))*((1/(2*h))*( theta(i+1-(N-1))-theta(i-1-(N-1)))).^3;

end
% for the two positons to the left of the rigbar ht hand boundary z=d
rhsode(2*N-3,1) = (B3(N-2)/(2*h^3))*( -theta(N-4) + 2*theta(N-3) - 2*theta(N-1) + Phi_b )                                + ((gbar (N-2)- C4(N-2))/h.^2)*(v(N-1) -2*v(N-2) + v(N-3))                  + ( A2(N-2)  +A3(N-2)  +2*B4(N-2) )*((theta(N-1)-theta(N-3))/(2*h)).*(theta(N-1) -2*theta(N-2) + theta(N-3))/(h^2)...
                  + ((Kappa1(N-2) -   2*C2(N-2))*(1/(2*h))*( v(N-1)-v(N-3) ) + D2(N-2)*sin(theta(N-2)).*cos(theta(N-2))  + E1(N-2)  + E2(N-2) )*(1/(2*h))*(theta(N-1)-theta(N-3))                + ( B2(N-2)  + C3(N-2) )*((1/(2*h))*(theta(N-1)-theta(N-3))).^3;

rhsode(2*N-2,1) = (B3(N-1)/(h^3)).*( -theta(N-3) + 3*theta(N-2) - 3*theta(N-1) + Phi_b )                                 + (( gbar (N-1)-C4(N-1))./(h.^2)).*( 0 -2*v(N-1) + v(N-2))                  + (A2(N-1) +A3(N-1) +2*B4(N-1) )* (( Phi_b-theta(N-2))/(2*h)).*( Phi_b -2*theta(N-1) + theta(N-2))/(h^2)...
                  + ( (Kappa1(N-1) - 2*C2(N-1))*(1/(2*h))*(0-v(N-2))   + D2(N-1).*sin(theta(N-1)).*cos(theta(N-1))       + E1(N-1) + E2(N-1))*(1/(2*h))*(  Phi_b-theta(N-2))                     + ( B2(N-1)  + C3(N-1) )*((1/(2*h))*( Phi_b-theta(N-2))).^3;
f = rhsode(N:2*N-2);
eqn = f == 0;
[A,b] = equationsToMatrix(eqn);
rank(A)
N-1
sol = A\b;
double(sol)
double(A*sol-b)
end 

[University]

Hi Torsten,

Thank you so much for your assistance. I will pick from here.