3D Plot for an integral function in a loop

Hello, I having some troubles plotting the function for the system below. Any assistance is greatly appreciated.

m1=5;
mtmd=0.25;
omg1=2*pi; %rad/s
beta=0.5
b=beta*m1
C1=0; 
K1=(omg1^2)*m1;

Freqratio=linspace(0.1,1,10);
zitad=linspace(0.1,1,10);
[zitad,Freqratio]=meshgrid(zitad,Freqratio);
out=zeros(size(zitad))

for f=0.1:0.1:length(Freqratio)
    for z=0.1:0.1:length(zitad)
     kTMDI=Freqratio(f)^2*omg1^2*(mtmd+b)
     cTMDI=2*zitad(z)*(mtmd+b)*Freqratio(f)*omg1

Ms=[mtmd 0;
    0 m1;]
M=Ms + [b 0;0 0;];
K=[kTMDI -kTMDI;
    -kTMDI kTMDI+K1;]
C=[cTMDI -cTMDI;
    -cTMDI cTMDI+C1;]

n=size(M,1)-1;
I=eye(2*n+2);
C0=zeros(2*n+2,1)' % Is (2n+2) length Vector
C0(2)=1

ae1=zeros(n+1);
ae2=eye(n+1);
ae3=-(M\K); %A\B= inv(B)*A;
ae4=-(M\C);
ae5=zeros(n+1,1);
ae6=ones(n+1,1);
A=[ae1 ae2;ae3 ae4];
B=[ae5;M\(Ms*ae6)];

whiteNoiseFactor=1;
g2=@(w) abs(C0*(((1i*w).*I-A)\B))^2*whiteNoiseFactor;
out(f, z)=integral(g2,0,inf,'ArrayValued',true);
    end 
end 

figure
surfc(Freqratio,zitad,out)

Answer 1

[VBBV]

Try now ... For loops cannot take decimal increments in matlab

%if true
  % code
%end
m1=5;
mtmd=0.25;
omg1=2*pi; %rad/s
beta=0.5
b=beta*m1
C1=0; 
K1=(omg1^2)*m1;
Freqratio=linspace(0.1,1,10);
zitad=linspace(0.1,1,10);
[zitad,Freqratio]=meshgrid(zitad,Freqratio);
out=zeros(size(zitad))
for f=1:length(Freqratio) % loops do not take decimal increments in matlab 
  for z=1:length(zitad)
  kTMDI=Freqratio(f)^2*omg1^2*(mtmd+b)
   cTMDI=2*zitad(z)*(mtmd+b)*Freqratio(f)*omg1
Ms=[mtmd 0;
  0 m1;]
M=Ms + [b 0;0 0;];
K=[kTMDI -kTMDI;
  -kTMDI kTMDI+K1;]
C=[cTMDI -cTMDI;
  -cTMDI cTMDI+C1;]
n=size(M,1)-1;
I=eye(2*n+2);
C0=zeros(2*n+2,1)' % Is (2n+2) length Vector
C0(2)=1
ae1=zeros(n+1);
ae2=eye(n+1);
ae3=-(M\K); %A\B= inv(B)*A;
ae4=-(M\C);
ae5=zeros(n+1,1);
ae6=ones(n+1,1);
A=[ae1 ae2;ae3 ae4];
B=[ae5;M\(Ms*ae6)];
whiteNoiseFactor=1;
g2=@(w) abs(C0*(((1i*w).*I-A)\B))^2*whiteNoiseFactor;
out(f, z)=integral(g2,0,inf,'ArrayValued',true);
  end 
end 
figure
surfc(Freqratio,zitad,out)% recommend to use mesh instead of surfc

[Komal Rajana]

Thank you. I have made the adjustment to the codes earlier but i am getting a different solution from the codes you modified. See below: I am not sure why is this?

m1=5;
mtmd=0.25;
omg1=2*pi; %rad/s
beta=0.5
b=beta*m1
C1=0; 
K1=(omg1^2)*m1;

F=0.1:0.1:1;
Fmax=numel(F)
z=0.1:0.1:1;
zmax=numel(z)

[z,F]=meshgrid(z,F);
out=zeros(size(F))

for i=1:Fmax
    for j=1:zmax
     kTMDI=F(i,j).^2*omg1^2.*(mtmd+b)
     cTMDI=2*z(i,j).*(mtmd+b).*F(i,j).*omg1

Ms=[mtmd 0;
    0 m1;]
M=Ms + [b 0;0 0;];
K=[kTMDI -kTMDI;
    -kTMDI kTMDI+K1;]
C=[cTMDI -cTMDI;
    -cTMDI cTMDI+C1;]

n=size(M,1)-1;
I=eye(2*n+2);
C0=zeros(2*n+2,1)' % Is (2n+2) length Vector
C0(2)=1

ae1=zeros(n+1);
ae2=eye(n+1);
ae3=-(M\K); %A\B= inv(B)*A;
ae4=-(M\C);
ae5=zeros(n+1,1);
ae6=ones(n+1,1);
A=[ae1 ae2;ae3 ae4];
B=[ae5;M\(Ms*ae6)];

whiteNoiseFactor=1;
g2=@(w) abs(C0*(((1i*w).*I-A)\B))^2*whiteNoiseFactor;
out(i,j)=integral(g2,-inf,inf,'ArrayValued',true);
    end 
end 
figure
surf(F,z,out)

[VBBV]

Because you changed the limits of integration in [0 inf] to [-inf inf]

 %if true
    % code
  %end
out(i,j)=integral(g2,-inf,inf,'ArrayValued',true)

You don't require a dot operator when using loop indices

%if true
  % code
% end
kTMDI=F(i,j)^2*omg1^2.*(mtmd+b)...

[Komal Rajana]

thanks very much.