plot function of two variables as one variable

I want to plot theta in the middle of the layer by setting z=d/2, for t=0..100. I have attached the snapshot of what I'm expecting. See my code below:

%% initialization
clear 
clc
close all

%% Model parameters
global k1 eta1 alpha3 gamma1 d N Phi xi h A B C G
k1 = 6*10^(-12);            %  elastic constant
eta1 = 0.0240;              %  viscosity
xi = 0.06;                  %  activity parameter
alpha3 = -0.001104;         %  viscosity
gamma1 = 0.1093;            %  viscosity
Theta = 0.0001;
% Simplify model parameters
A = k1/gamma1;
B = alpha3/gamma1;
C = eta1 - alpha3*B;
G = alpha3*A;
d = 0.0002;
N = 20;

%% Numerical setup
% step size
h = d/N;
% t span
tspan = [0 100];
% range of z
%z=linspace(0,d/2,N+1)
z=linspace(0,d,N+1);
% initial conditions
theta0 = Theta*sin(pi*z/d);
v0 = zeros(1,N+1);
theta0_int=theta0(2:N);
v0_int=v0(2:N);
u0 = [theta0_int'; v0_int'];
% Constant mass matrix M: We use M to represent the left hand side of the system of equations.
M1=eye(N-1,N-1);
M2=zeros(N-1,N-1);
M=[M1 M2;M2 M2];
% Boundary Conditions
Phi = 0;
%% ode solver
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
[t,y] = ode15s(@lcode1,tspan,u0, options);
% Extract the solution for theta and v
theta = [Phi*ones(length(t), 1) y(:,1:N-1) Phi*ones(length(t), 1)];
v = [zeros(length(t), 1) y(:,N:(2*N-2)) zeros(length(t), 1)];
%% Plot the solution.
figure
plot(z, theta(end,:))
xlabel('z(\mum)')
ylabel('\theta(z,t)(rad)')
title('Director angle plot')
%% Functions
% Right hand side of the ODEs: F(U)
function rhsode = lcode1(t,y)
global Phi xi h A B C G N
% initialize theta and v
theta = y(1:(N-1));
v = y(N:(2*N-2));
% theta equations
% theta equations
% for the positon to the right of the left hand boundary z=0
rhsode(1,1) =     (A/(h^2))*(theta(2)-2*theta(1)+ Phi)          - (B/(2*h))*(v(2)-0);
% for all other internal positions 0<z<d
for i=2:(N-2)
    rhsode(i,1) = (A/(h^2))*(theta(i+1)-2*theta(i)+theta(i-1))  - (B/(2*h))*(v(i+1)-v(i-1)); 
end 
% for the positon to the left of the right hand boundary z=d
rhsode(N-1,1)   = (A/(h^2))*(Phi-2*theta(N-1)+ theta(N-2))      - (B/(2*h))*(0-v(N-2));
% 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 right of the left hand boundary z=0
rhsode(N,1)     = (G/(h^3))*(-Phi + 3*theta(1) - 3*theta(2) + theta(3))            + (C/(h^2))*(v(2) -2*v(1)+ 0)             + (xi/(2*h))*(theta(2)-Phi);
rhsode(N+1,1)   = (G/(2*h^3))*(Phi - 2*theta(1) + 2*theta(3)- theta(4))            + (C/(h^2))*(v(3) -2*v(2) + v(1))         + (xi/(2*h))*(theta(3)-theta(1));
% for all other internal positions 0<z<d
for i=(N+2):(2*N-4)
    rhsode(i,1) = (G/(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)))   +(C/(h^2))*(v(i+1-(N-1)) -2*v(i-(N-1)) + v(i-1-(N-1)))          + (xi/(2*h))*(theta(i+1-(N-1))-theta(i-1-(N-1)));
end
% for the two positons to the left of the right hand boundary z=d
rhsode(2*N-3,1) = (G/(2*h^3))*(theta(N-4) - 2*theta(N-3) + 2*theta(N-1) - Phi)     +(C/(h^2))*(v(N-1) -2*v(N-2) + v(N-3))    + (xi/(2*h))*(theta(N-1)-theta(N-3));
rhsode(2*N-2,1) = (G/(h^3))*(-theta(N-3) + 3*theta(N-2) - 3*theta(N-1) + Phi)      +(C/(h^2))*(0 -2*v(N-1) + v(N-2))         + (xi/(2*h))*(Phi-theta(N-2));      
end

Answer 1

[Torsten]

%% initialization
clear 
clc
close all
%% Model parameters
global k1 eta1 alpha3 gamma1 d N Phi xi h A B C G
k1 = 6*10^(-12);            %  elastic constant
eta1 = 0.0240;              %  viscosity
xi = 0.06;                  %  activity parameter
alpha3 = -0.001104;         %  viscosity
gamma1 = 0.1093;            %  viscosity
Theta = 0.0001;
% Simplify model parameters
A = k1/gamma1;
B = alpha3/gamma1;
C = eta1 - alpha3*B;
G = alpha3*A;
d = 0.0002;
N = 20;
%% Numerical setup
% step size
h = d/N;
% t span
tspan = [0 230];
% range of z
%z=linspace(0,d/2,N+1)
z=linspace(0,d,N+1);
% initial conditions
theta0 = Theta*sin(pi*z/d);
v0 = zeros(1,N+1);
theta0_int=theta0(2:N);
v0_int=v0(2:N);
u0 = [theta0_int'; v0_int'];
% Constant mass matrix M: We use M to represent the left hand side of the system of equations.
M1=eye(N-1,N-1);
M2=zeros(N-1,N-1);
M=[M1 M2;M2 M2];
% Boundary Conditions
Phi = 0;
%% ode solver
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
[t,y] = ode15s(@lcode1,tspan,u0, options);
theta_middle = y(:,N/2);
plot(t,theta_middle)

% Extract the solution for theta and v
%theta = [Phi*ones(length(t), 1) y(:,1:N-1) Phi*ones(length(t), 1)];
%v = [zeros(length(t), 1) y(:,N:(2*N-2)) zeros(length(t), 1)];
%% Plot the solution.
%figure
%plot(z, theta(end,:))
%xlabel('z(\mum)')
%ylabel('\theta(z,t)(rad)')
%title('Director angle plot')
%% Functions
% Right hand side of the ODEs: F(U)
function rhsode = lcode1(t,y)
global Phi xi h A B C G N
% initialize theta and v
theta = y(1:(N-1));
v = y(N:(2*N-2));
% theta equations
% theta equations
% for the positon to the right of the left hand boundary z=0
rhsode(1,1) =     (A/(h^2))*(theta(2)-2*theta(1)+ Phi)          - (B/(2*h))*(v(2)-0);
% for all other internal positions 0<z<d
for i=2:(N-2)
    rhsode(i,1) = (A/(h^2))*(theta(i+1)-2*theta(i)+theta(i-1))  - (B/(2*h))*(v(i+1)-v(i-1)); 
end 
% for the positon to the left of the right hand boundary z=d
rhsode(N-1,1)   = (A/(h^2))*(Phi-2*theta(N-1)+ theta(N-2))      - (B/(2*h))*(0-v(N-2));
% 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 right of the left hand boundary z=0
rhsode(N,1)     = (G/(h^3))*(-Phi + 3*theta(1) - 3*theta(2) + theta(3))            + (C/(h^2))*(v(2) -2*v(1)+ 0)             + (xi/(2*h))*(theta(2)-Phi);
rhsode(N+1,1)   = (G/(2*h^3))*(Phi - 2*theta(1) + 2*theta(3)- theta(4))            + (C/(h^2))*(v(3) -2*v(2) + v(1))         + (xi/(2*h))*(theta(3)-theta(1));
% for all other internal positions 0<z<d
for i=(N+2):(2*N-4)
    rhsode(i,1) = (G/(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)))   +(C/(h^2))*(v(i+1-(N-1)) -2*v(i-(N-1)) + v(i-1-(N-1)))          + (xi/(2*h))*(theta(i+1-(N-1))-theta(i-1-(N-1)));
end
% for the two positons to the left of the right hand boundary z=d
rhsode(2*N-3,1) = (G/(2*h^3))*(theta(N-4) - 2*theta(N-3) + 2*theta(N-1) - Phi)     +(C/(h^2))*(v(N-1) -2*v(N-2) + v(N-3))    + (xi/(2*h))*(theta(N-1)-theta(N-3));
rhsode(2*N-2,1) = (G/(h^3))*(-theta(N-3) + 3*theta(N-2) - 3*theta(N-1) + Phi)      +(C/(h^2))*(0 -2*v(N-1) + v(N-2))         + (xi/(2*h))*(Phi-theta(N-2));      
end

[University Glasgow]

Thank you very much.

[University Glasgow]

Please how do I save figure without white space?

[Torsten]

https://de.mathworks.com/help/matlab/creating_plots/save-figure-with-minimal-white-space.html

[University Glasgow]

I tried this ax = gca;

% Requires R2020a or later

exportgraphics(ax,'myplot.png','Resolution',300) but is saving only one of my plot.

%% initialization

clear

clc

close all

%% Model parameters

global k1 eta1 alpha3 gamma1 d N Phi xi h A B C G

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

eta1 = 0.0240; % viscosity

xi = -1; % activity parameter

alpha3 = -0.001104; % viscosity

gamma1 = 0.1093; % viscosity

Theta = 0.0001;

% Simplify model parameters

A = k1/gamma1;

B = alpha3/gamma1;

C = eta1 - alpha3*B;

G = alpha3*A;

d = 0.0002;

N = 200;

%% Numerical setup

% step size

h = d/N;

% t span

tspan = [0 100];

% range of z

%z=linspace(0,d/2,N+1)

z=linspace(0,d,N+1);

% initial conditions

theta0 = Theta*sin(pi*z/d);

v0 = zeros(1,N+1);

theta0_int=theta0(2:N);

v0_int=v0(2:N);

u0 = [theta0_int'; v0_int'];

% Constant mass matrix M: We use M to represent the left hand side of the system of equations.

M1=eye(N-1,N-1);

M2=zeros(N-1,N-1);

M=[M1 M2;M2 M2];

% Boundary Conditions

Phi = 0;

%% ode solver

options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);

[t,y] = ode15s(@lcode1,tspan,u0, options);

theta_middle = y(:,N/2); % theta at the middle of the layer (i.e., z=d/2)

% Extract the solution for theta and v

theta = [Phi*ones(length(t), 1) y(:,1:N-1) Phi*ones(length(t), 1)];

v = [zeros(length(t), 1) y(:,N:(2*N-2)) zeros(length(t), 1)];

%% Plot the solution.

figure

subplot(2,2,1)

plot(z, theta(1:5,:))

xlabel('z(\mum)')

ylabel('\theta(z,t)(rad)')

title('Director angle','FontSize',10)

%[xmin(z) xmax(z) ymin(theta) ymax(theta)];

subplot(2,2,2)

plot(z, v(1:5,:))

xlabel('z(\mum)')

ylabel('v(z,t)(m/s)')

title('Flow velocity')

%[xmin(z) xmax(z) ymin(theta) ymax(theta)];

subplot(2,2,3)

plot(t,theta_middle)

xlabel('z(\mum)')

ylabel('\theta(d/2,t)(rad)')

title('Director angle at the middle of the layer', 'FontSize',6)

ax = gca;

% Requires R2020a or later

exportgraphics(ax,'modelplot.png','Resolution',300)

%% Functions

% Right hand side of the ODEs: F(U)

function rhsode = lcode1(t,y)

global Phi xi h A B C G N

% initialize theta and v

theta = y(1:(N-1));

v = y(N:(2*N-2));

rhsode(1,1) = (A/(h^2))*(theta(2)-2*theta(1)+ Phi) - (B/(2*h))*(v(2)-0);

% for all other internal positions 0<z<d

for i=2:(N-2)

rhsode(i,1) = (A/(h^2))*(theta(i+1)-2*theta(i)+theta(i-1)) - (B/(2*h))*(v(i+1)-v(i-1));

end

% for the positon to the left of the right hand boundary z=d

rhsode(N-1,1) = (A/(h^2))*(Phi-2*theta(N-1)+ theta(N-2)) - (B/(2*h))*(0-v(N-2));

rhsode(N,1) = (G/(h^3))*(-Phi + 3*theta(1) - 3*theta(2) + theta(3)) + (C/(h^2))*(v(2) -2*v(1)+ 0) + (xi/(2*h))*(theta(2)-Phi);

rhsode(N+1,1) = (G/(2*h^3))*(Phi - 2*theta(1) + 2*theta(3)- theta(4)) + (C/(h^2))*(v(3) -2*v(2) + v(1)) + (xi/(2*h))*(theta(3)-theta(1));

% for all other internal positions 0<z<d

for i=(N+2):(2*N-4)

rhsode(i,1) = (G/(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))) +(C/(h^2))*(v(i+1-(N-1)) -2*v(i-(N-1)) + v(i-1-(N-1))) + (xi/(2*h))*(theta(i+1-(N-1))-theta(i-1-(N-1)));

end

% for the two positons to the left of the right hand boundary z=d

rhsode(2*N-3,1) = (G/(2*h^3))*(theta(N-4) - 2*theta(N-3) + 2*theta(N-1) - Phi) +(C/(h^2))*(v(N-1) -2*v(N-2) + v(N-3)) + (xi/(2*h))*(theta(N-1)-theta(N-3));

rhsode(2*N-2,1) = (G/(h^3))*(-theta(N-3) + 3*theta(N-2) - 3*theta(N-1) + Phi) +(C/(h^2))*(0 -2*v(N-1) + v(N-2)) + (xi/(2*h))*(Phi-theta(N-2));

end

[Torsten]

Sorry, but I have no experience in optimizing graphics output.

Maybe you should open a new question.