Need help plotting this

I need help plotting this code so that my plot will look like the discontinuous sine waves (bolded black) in the graph posted below the code (for one unit cell). Right now my plot is showing up blank.

clear all;
format long;
im = sqrt(-1);
CellLength = 1;
ibeta = 1;

%Define materal properties 
CellLength = 1;
layers = 2;
d = [0.4;0.6];
dTotal = d(1,1)+d(2,1);

xc = [0;0.4];

Ef = 12;
pf = 3;
cf = sqrt(Ef/pf);

Em = 1;
pm = 1;
cm = sqrt(Em/pm);

w = 5;

T1 = [cos(.2*w) (1/(6*w))*sin(.2*w); -6*w*sin(.2*w) cos(.2*w)];
T2 = [cos(.6*w) (1/w)*sin(.6*w); -w*sin(.6*w) cos(.6*w)];
T = T2*T1;

Z1 = 6*w;
Z2 = w;
Z = [Z1;Z2];

%Solve eigenvalue problem for k
[V,D] = eig(T);                 %D = eigenvalues, %V = eigenvectors
k1 = log(D(1,1))/(im*dTotal);
k2 = log(D(2,2))/(im*dTotal);
k = [k1;k2];

B1 = [1 1;im*30 -im*30];
B2 = [1 1;im*5 -im*5];
C1a = [1 0;0 1];
C2a = [exp(im*k2*0.4) 0;0 exp(-im*k2*0.4)];
a1 = inv(B1)*V(:,1);
a2 = inv(C2a)*inv(B2)*T1*B1*a1;

for x1 = 0:0.1:0.4
    C1 = @(x1)([exp(im*k1*x1) 0;0 exp(-im*k1*x1)]);
    C1 = C1(x1);
    y1 = @(x1)(B1*C1*a1);
    y1 = y1(x1);
end
for x2 = 0.4:0.1:1
    C2 = @(x2)([exp(im*k2*x2) 0;0 exp(-im*k2*x2)]);
    C2 = C2(x2);
    y2 = @(x2)(B2*C2*a2);
    y2 = y2(x2);
end

plot(x1,real(y1(1,:)))
hold on
plot(x2,real(y2(1,:)))
xlabel('Position, x')
ylabel('Displacement, u')

Answer 1

[Star Strider]

Two related problems are that the ‘y1’ and ‘y2’ functions do not use their arguments in their calculations. (I would also rename them ‘y1f’ and ‘y2f’ to avoid confusion with your ‘y1’ and ‘y2’ vectors.)

What do you want to do in those functions?

[Amanda Lococo]

I wasn't sure if I needed the functions for those or not. I figured I didn't since the only x-iteration is in the function C2. I am trying to write a code based on the results of a paper and the function described is that which is stated: y(x) = BjCj(x)aj to be specific.

What I need the function to do overall is plot x against y, where y is changing as a function of x (because C depends on x).

[Star Strider]

I have no idea what you’re doing. This is a bit more efficient, and produces plots:

ibeta = 1;
%Define materal properties 
CellLength = 1;
layers = 2;
d = [0.4;0.6];
dTotal = d(1,1)+d(2,1);

xc = [0;0.4];

Ef = 12;
pf = 3;
cf = sqrt(Ef/pf);

Em = 1;
pm = 1;
cm = sqrt(Em/pm);

w = 5;

T1 = [cos(.2*w) (1/(6*w))*sin(.2*w); -6*w*sin(.2*w) cos(.2*w)];
T2 = [cos(.6*w) (1/w)*sin(.6*w); -w*sin(.6*w) cos(.6*w)];
T = T2*T1;

Z1 = 6*w;
Z2 = w;
Z = [Z1;Z2];

%Solve eigenvalue problem for k
[V,D] = eig(T);                 %D = eigenvalues, %V = eigenvectors
k1 = log(D(1,1))/(1i*dTotal);
k2 = log(D(2,2))/(1i*dTotal);
k = [k1;k2];

B1 = [1 1;1i*30 -1i*30];
B2 = [1 1;1i*5 -1i*5];
C1a = [1 0;0 1];
C2a = [exp(1i*k2*0.4) 0;0 exp(-1i*k2*0.4)];
a1 = B1\V(:,1);
a2 = (B2*C2a)\T1*B1*a1;

x1 = 0:0.1:0.4;
C1 = @(x1)([exp(1i*k1*x1) 0;0 exp(-1i*k1*x1)]);
y1 = zeros(2, numel(x1));
for k = 1:numel(x1)
    y1(:,k) = (B1*C1(x1(k))*a1);
end

x2 = 0.4:0.1:1;
C2 = @(x2)([exp(1i*k2*x2) 0;0 exp(-1i*k2*x2)]);
y2 = zeros(2, numel(x2));
for k = 1:numel(x2)
    y2(:,k) = (B2*C2(x2(k))*a2);
end

plot(x1,real(y1(1,:)))
hold on
plot(x2,real(y2(1,:)))
xlabel('Position, x')
ylabel('Displacement, u')

I can’t completely vectorise this because I’m not certain what you’re doing. I moved the function definitions out of the loops, and called them in the ‘y1’ and ‘y2’ calculations.

You’ll have to determine if the plots are correct. I added preallocations for ‘y1’ and ‘y2’.

[Amanda Lococo]

The plots do work, but the plots I need are posted in the original thread. It should resemble a sinusoidal wave with discontinuities at the layer boundaries (in this case at 0.4). I appreciate all of your help!

[Star Strider]

My pleasure.

If you extend ‘x2’ to at least 3, you get the sine curve:

x2 = 0.4:0.1:3;

and with ‘x1’ defined as it is, the discontinuity appears at 0.4.

Now that your code produces plots, you need to experiment with it to get the results you want. I have no idea what you are doing, so I can only help with the code, not the concept (that may be outside my areas of expertise).

[Amanda Lococo]

I will try that then! Thanks again for all of your help!

[Star Strider]

As always, my pleasure!