function sol_nonl_eq
clear all;clc
global x
format long
y0=0.1;
xspan=0.1:0.1:1;
n=length(xspan);
for i=1:n
    x=xspan(i);
    y(i)=fsolve(@fun,y0);
    y_fval(i)=fun(y(i));
end
[xspan' y' y_fval']

plot(xspan,y,'b-o')

   


function f=fun(y)
global x
f=y-3.5e-10/(x^2)/(10^(-2.0588*(x+3*y)^0.5/(1+1.9758*(x+3*y)^0.5)))/(10^(-0.5147*(x+3*y)^0.5/(1+4.1393*(x+3*y)^0.5)))^2;