Using fsolve after matlabFunction

Hello everyone! Hope you are doing well.

I was getting stuck at a point while coding in MATLAB. I have to solve an equation which looks like as shown below -

function [vr,vt]=Stagpt(r,theta,Q,N_t) % Q-flow rate 

n=1:N_t;
m=5;
An=0;
Bn=0;
for j=1:m
    An=An+(2./(n*pi)).*(Q(j)*sin(n*(6-j)*pi/6).*cos(n*(4+j)*pi/6));
    Bn=Bn+(2./(n*pi)).*(Q(j)*sin(n*(6-j)*pi/6).*sin(n*(4+m)*pi/6));
end
vr=0;
vt=0;
for i=1:N_t
    vr=vr+((i*r^(i-1)).*(-An(i)*sin(i*theta)+Bn(i)*cos(i*theta))); % In both these expressions r and theta are symbolic variables.
    vt=vt-((i*r^(i-1)).*(An(i)*sin(i*theta)+Bn(i)*cos(i*theta)));  
end

Now the main function for my case looks like as follows

syms r theta
[vr,vt]= Stagpt(r,theta,[1 -1 1 -1 1],10);

v = matlabFunction(vr,vt,'File','Velocity','Optimize',false);

After giving the matlabFunction command I get the following function

function [vr,vt] = Velocity(r,theta)
%VELOCITY
%    [VR,VT] = VELOCITY(R,THETA)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    26-Jan-2019 13:33:43

vr = cos(theta).*(-1.705817538915785e-1)-sin(theta).*5.254972497472503e-17+r.^5.*(cos(theta.*6.0).*4.225610904505511e-31+sin(theta.*6.0).*2.27459032782291e-16).*6.0+(r.*sin(theta.*2.0))./2.251799813685248e15-r.^6.*(cos(theta.*7.0).*3.394139050826781e-1+sin(theta.*7.0).*1.333257168849096e-16).*7.0+r.^8.*(cos(theta.*9.0).*7.073553026306457e-2-sin(theta.*9.0).*6.925704865457059e-17).*9.0+r.^3.*(cos(theta.*4.0).*1.109335647967048e-31+sin(theta.*4.0).*2.775557561562891e-16).*4.0+r.^7.*(cos(theta.*8.0).*2.218671295934096e-31+sin(theta.*8.0).*2.359223927328458e-16).*8.0-r.^9.*(cos(theta.*1.0e1).*3.57452597678271e-31+sin(theta.*1.0e1).*7.632783294297951e-17).*1.0e1+r.^2.*(cos(theta.*3.0).*2.122065907891938e-1-sin(theta.*3.0).*3.83385906083066e-16).*3.0-r.^4.*(cos(theta.*5.0).*4.751794671157494e-1+sin(theta.*5.0).*4.76749172825605e-17).*5.0;
if nargout > 1
    vt = cos(theta).*1.705817538915785e-1-sin(theta).*5.254972497472503e-17-r.^5.*(cos(theta.*6.0).*4.225610904505511e-31-sin(theta.*6.0).*2.27459032782291e-16).*6.0+(r.*sin(theta.*2.0))./2.251799813685248e15+r.^6.*(cos(theta.*7.0).*3.394139050826781e-1-sin(theta.*7.0).*1.333257168849096e-16).*7.0-r.^8.*(cos(theta.*9.0).*7.073553026306457e-2+sin(theta.*9.0).*6.925704865457059e-17).*9.0-r.^3.*(cos(theta.*4.0).*1.109335647967048e-31-sin(theta.*4.0).*2.775557561562891e-16).*4.0-r.^7.*(cos(theta.*8.0).*2.218671295934096e-31-sin(theta.*8.0).*2.359223927328458e-16).*8.0+r.^9.*(cos(theta.*1.0e1).*3.57452597678271e-31-sin(theta.*1.0e1).*7.632783294297951e-17).*1.0e1-r.^2.*(cos(theta.*3.0).*2.122065907891938e-1+sin(theta.*3.0).*3.83385906083066e-16).*3.0+r.^4.*(cos(theta.*5.0).*4.751794671157494e-1-sin(theta.*5.0).*4.76749172825605e-17).*5.0;
end

Now, my question is I want to simultaneously solve the expressions vr==0 and vt==0 so that I get r and theta values respectively. One solution is to copy past the

expression and use fsolve, but I also want to change the value of Q later (give it as an input by the user). I know I have to use fsolve but I am not sure how to proceed.

Can someone please help me with this? Or just direct me to the correct command or another question of the same kind.

Thanking you!