The problem is about initial guess of your equation. vpasolve can not correctly find an initial point and therefore returns a wrong value. My suggestion would be take a look of the graph of the function. First, run the following code:
syms yn
b=3;
m=1;
S=0.001;
n=0.014;
Q=28;
A=(b+m*yn)*yn;
T=b+2*m*yn;
P=b+2*yn*sqrt(1+m^2);
R=A/P;
D=A/T;
fun=Q-(A./n).*S.^(0.5).*R.^(2/3);
plot(-20:0.1:20,subs(fun,yn,-20:0.1:20),-20:0.1:20,subs(fun,yn,-20:0.1:20)==0,'*')
When you check the plot for the input vector and when you zoom in the red line crosses the function, you will see that it crosses when yn has a value of 2.1 or around it. Therefore, you need to set your initial guess argument accordingly, so that you will have an accurate solution. After this step, run the following code and see the result.
yn=vpasolve(fun,yn,3)
I chose 3 as initial point in this situation.
