Why can't I get matlab to do this when mathematica can do it easily

I have written code in mathematica which can successfully generate a figure for me, and I'm trying to do the same thing in matlab but I've been stuck on it for a while now. I'm pretty sure it is my understanding of how matlab operates that is the problem, (as opposed to matlab's inability to do what mathematica can do) but I haven't been able to make it work.

Could someone take a look at my code and determine what I'm doing wrong with the code, it is actually pretty simple. Look at the matematica code first to see what it should be doing.

I've attached 4 pictures (mathematica 1 & 2 ....also Matlab 1 & 2): Two are of the mathematica code/plot. The first mathematica pic being the first function plotted and then the second mathematica pic is the end result. And there are two pics of my matlab code/plot - the first is of the first function (which is identical to the mathematica one according to the graph), and then the second pic is the end result (which is wrong when plotting EOM2 according to the graph).

----------------------------------------------------------------------------------------------

My mathematica code looks like this:

f[x_] := PDF[NormalDistribution[], x]
v1 := 1.1
v2 := .9
f1 := 13.3
f2 := 20
EOM1[x_] := Sum[BesselJ[Abs[i], v1]*f[x + f1*i], {i, -2, 2}]
EOM2[x_] := Sum[BesselJ[Abs[i], v2]*EOM1[x + f2*i], {i, -2, 2}]

Overlay[{ 
  Plot[EOM1[x], {x, -100, 100}, PlotRange -> {0, .5}, 
   PlotStyle -> Red],
  Plot[12*PDF[NormalDistribution[0, 20], x], {x, -100, 100},
   PlotRange -> {0, .5}]
  }]

-----------------------------------------------------------------------------------------------

My MATLAB code looks like this:

clear all; clc;
xStep = 1;
x = (-50:xStep:50);
V1 = 1.1;
V2 = .9;
F1 = 13.3;
F2 = 20;
%%First part

k = (-2:1:2);
s = 1;
u = 0;
for x1 = min(x):max(x)
    for  k1 = min(k):max(k)

          %BSL1(k1+max(k)+1,x1+max(x)+1) = besselj(abs(k1),V1).*normpdf((x1+F1.*k1),0,1);
          BSL2a(k1+max(k)+1,x1+max(x)+1) = besselj(abs(k1),V1).*(1/(s*sqrt(2*pi))).*exp((-1./2*s.^2).*(x1+F1.*k1-u).^2);

      end
  end
  BSL2aSum = sum(BSL2a);
  plot(x,BSL2aSum,'b'); hold off;

%%Second part...
for x1 = min(x):max(x)
    for  k1 = min(k):max(k)

          %EOM2(k2+max(k)+1,k1+max(k)+1,x1+max(x)+1) = besselj(abs(k1),V1).*besselj(abs(k2),V2).* (1/(s*sqrt(2*pi))).*exp((-1./2*s.^2).*(x1+F1.*k1+F2.*k2-u).^2);
          EOM2(k1+max(k)+1,x1+max(x)+1)  = normpdf((x1+(F1.*k1)+(F2.*k1)),0,1).*besselj(abs(k1),V1).*besselj(abs(k1),V2);

      end
  end
  plot(x,sum(EOM2))

Answer 1

[caleb]

That's it, you nailed it Joesph! Thanks!

As you can see from one of the lines commented out in my matlab code, I had tried that before, but the order of my 3 for loops was wrong. I re-ordered them so that x is looped through first then the 2 other k-variables and that fixed my problem.

Answer 2

[Joseph Cheng]

What is happening is you're not getting the double loop through k1 for EOM2. If you look at the Mathematica code:

EOM1[x_] := Sum[BesselJ[Abs[i], v1]*f[x + f1*i], {i, -2, 2}]
EOM2[x_] := Sum[BesselJ[Abs[i], v2]*EOM1[x + f2*i], {i, -2, 2}]

EOM2(y) has EOM1(y+f2*i) for i=-2 to +2.

Well if we go one step deeper then f(x+f1*i) for i=-2 to +2 in EOM1. But here the x in the f(x) is actually replaced with the y from -2 to 2. hopefully that make sense... i'll re-write it:

in EOM2 y+f2*n for n=-2:2 and then working it into EOM1 y+f2n+f1*i for n=-2:2 for i=-2:2. (then you see you're missing the n=-2 for all cases of i, n=-1 for all cases of i, etc)

[Joseph Cheng]

so you got very close...

for x1 = min(x):max(x)
    for  k1 = min(k):max(k)
        for k2 = min(k):max(k)
          EOM2(k2+max(k)+1,k1+max(k)+1,x1+max(x)+1) = besselj(abs(k1),V1).*besselj(abs(k2),V2).* (1/(s*sqrt(2*pi))).*exp((-1./2*s.^2).*(x1+F1.*k1+F2.*k2-u).^2);
        end
      end
  end
  plot(x,squeeze(sum(sum(EOM2))))

i just saw your commented out code and the first attempt was actually good. but to plot it you needed to do the sum(sum(EOM2)) as you did get the 5x5xN where you needed to do the sum twice and then make it a 1xN for plot.

[Joseph Cheng]

Additionally if you want something that looks a bit more like the mathematica you can us the anonymous functions to define the two functions.

myfun1 = @(x,f1,i,V1,s,u) besselj(abs(i),V1).*(1/(s*sqrt(2*pi))).*exp((-1./2*s.^2).*(x+F1.*i-u).^2);
myfun2 = @(x,f1,f2,i,V1,V2,s,u)  besselj(abs(i),V2).*myfun1(x+f2.*i,f1,i,V1,s,u);

for x1=1:length(x)
    for k1=1:length(k)
        X1(x1,k1) = myfun1(x(x1),F1,k(k1),V1,s,u);
    end
end
X1 =sum(X1,2);
for x1=1:length(x)
    for k1=1:length(k)
        for k2=1:length(k)
            X2(k1,k2,x1) = besselj(abs(k(k2)),V2).*myfun1(x(x1)+F2.*k(k2),F1,k(k1),V1,s,u);
        end
    end
end
X2 = squeeze(sum(sum(X2)));

figure,subplot(2,1,1),plot(x,X1)
subplot(2,1,2),plot(x,X2)

[caleb]

you're the man! thanks for you help... I meant to select your first response as the correct answer but somehow mine got clicked. If "best answers" mean anything to you and you know how to correct that, please let me know and I'll switch it to your's.