I think what you are trying to do is vary Vds and see what Ids becomes. So here is a simple some modification for it using for loop. However you should try to vectorize the for loop for efficient code:
Vth=2e-7;
Vs=1.5e-7;
Z=150e-6;
e=1.6e-19;
ep=8.85e-12;
mn=9.1e-31;
dd=30e-9;
di=7e-9;
x=1.05e-34;
u=1.257e-6;
L=150e-6;
Vgs=-4;
Vds=-60:0.01:60; % generate a vector of interest
Vdsat=Vgs-Vth+Vs-((Vgs-Vth)*(Vgs-Vth)+Vs*Vs);
for i = 1:length(Vds) % loop through Vds and calculate corresponding Ids
if(Vdsat<Vds(i))
Ids(i)=2*Z*e*e*ep*mn*Vs/(2*(dd+di)*e*e*mn+ep*pi*x*x)*(((Vgs-Vth)*(Vgs-Vth)+Vs^0.5)-Vs);
end
if(Vdsat>Vds(i))
Ids(i)=2*u*Z*ep*e*e*mn/(L*(2*(dd+di)*e*e*mn+ep*pi*x*x))*((Vgs-Vth)*Vds(i)-((Vds(i).*Vds(i))/2));
end
end
figure;
plot(Vds,Ids);
xlabel('Vds (V)');
ylabel('Ids (mA/mm)');
