There are several issues, which all need to be settled:

You differentiate w.r.t. x, but no x appears anywhere.
What is the unknown function? One should think it would be F(s), but F(s) is defined to be a specific expression in s.
In Maple you cannot use square brackets in lieu of parentheses.
The notation for muo and muw change through the worksheet.

There may be more.

U as such is not known ever, but V is defined to be the value of U at tk.
Change 'if U(tk) <=4 to 'V <= 4 ´.

restart;
k := 0; B[1] := 0.001/365; B[2] :=0.002/365:
ode := diff(U(t), t) = -(A[1]+A[2]*U(t))*U(t);
ic[0] := U(365*k) = 1000;
sol[0] := dsolve({ic[0], subs(A[1] = B[1], A[2] = B[2], ode)}, U(t), numeric);
sigma := .5;
for k to 4 do
tk := 365*k;
V := rhs(sol[k-1](tk)[2]);
if  V <= 4 then ic[k] := U(tk) =0
else ic[k] := U(tk) = sigma*V
end if;
sol[k] := dsolve({ic[k], subs(A[1] = B[1], A[2] = B[2], ode)}, U(t), numeric)
end do;

#Plot:

for k to 5 do p[k-1]:=plots:-odeplot(sol[k-1],[t,U(t)],365*(k-1)..365*k) end do:
plots:-display(seq(p[k],k=0..4),thickness=2);

To see the effect of the 'if' clause you can replace the 4 by e.g. 40.

@argand Since this is an indefinite integral you have to allow for a constant. In this case it is imaginary!

int(U, t);
evalc(%);
combine(%);
evalc(%);
evalc(Im(%));
normal(%);
#Thus the imaginary part is constant!

Using evalf helps:

LinearAlgebra:-Eigenvalues(J(25/8, 5/18, -5/8));
evalf(%);
You will see that 2 of the eigenvalues are imaginary and have positive real parts. The third is (real and) negative.

For a bunch of orbits you could try something like this:

interface(warnlevel=0); #Initially you may not want to turn warnings off as they do give you some valuable information.
bg:=pointplot3d(pts,symbol=solidsphere,symbolsize=20,color=blue):
x0:=-2: y0:=-1: z0:=-2: dx:=6: dy:=3: dz:=3: N:=4:
t1:=-10: t2:=5:
DEplot3d({sys},[x(t),y(t),z(t)],t=t1..t2,[seq(seq(seq([x(0)=x0+i*dx/N,y(0)=y0+j*dy/N,z(0)=z0+k*dz/N],k=1..N-1),j=1..N-1),i=1..N-1)],linecolor=red,x=x0..x0+dx,y=y0..y0+dy,z=z0..z0+dz,thickness=2,stepsize=.01):
display(%,bg);

Admittedly a rather confusing picture.

When I pasted your procedure into a worksheet there was a multiplication sign following immediately after 'irem':

if evalb(irem*(u, 2) = 0)

That `*` certainly shouldn't be there.

In the Windows version of Maple the classical interface has always been included (also in Maple 16). I don't know about Macs.

However, the standard interface in (1) worksheet mode and (2) with Maple input (1D input) is (I think) fine.
You can make those 2 changes via the menu. Since I don't have a Mac I cannot go into the details. In Windows it is via Tools/Options, but I think it may be different on a Mac.

The following process finds a solution satisfying the initial condition, but it does not satisfy the boundary condition.

restart;
pde:=diff(u(x,t),t)+diff(u(x,t),x)+u(x,t)=Dirac(x);
ans:=pdsolve(pde);
eval(rhs(ans),t=0)=100;
solve(%,{_F1(-x)});
subs(x=x-t,%);
eval(ans,%);
res:=combine(expand(%));
#Testing the solution (a zero result means that it passes the test):
pdetest(res,pde);
eval(res,t=0);
eval(res,x=0);
#The 'undefined' part is due to H(0) not being defined in Maple. This can be seen by using subs instead:
subs(x=0,res);

(1/3)*``(3*a);

expand knows ``(), so expand gets you back to where you were:

expand(%);

I used a slightly different approach. SSQ is the sum of squares of the differences between the iV values and the i-values from the model. You will see that your initial guess is not very good.

(That may account for the problems NonlinearFit is having, but maybe not quite. NonlinearFit doesn't even return a result (in the time I allotted) when given the parameter values found below as initialvalues (with maxfun=0 in dsolve). Without maxfun=0  NonlinearFit complained that it couldn't go beyond 12.276 and the parameters obtained at that time were beta = 0.028879..., alpha = -.59691..., which are definitely lousy.)
Please also see addendum at bottom!

With everything above N defined as before, you can do as follows:

N := dsolve(sys, numeric, parameters = [beta, alpha], output = listprocedure);
ip:=subs(N,i(t));
SSQ:=proc(beta,alpha) local j;
     N(parameters=[beta,alpha]);
     add((ip(j)-iV[j])^2,j=1..len)
end proc;
SSQ(.5,.3);
plot3d(SSQ,0..0.2,0..0.5,axes=boxed,labels=[beta,alpha,SSQ]);
plot3d(SSQ,0..0.01,0..0.5,axes=boxed,labels=[beta,alpha,SSQ]);
contourplot(SSQ,0.002..0.008,0.2..0.5,labels=[beta,alpha],contours=50);
res:=Minimize(SSQ,0.002..0.008,0.2..0.5);
res[2];
N(parameters=convert(res[2],list));

odeplot(N,[t,i(t)],1..len,color=blue):
display(%,pp);

Addendum.
Again with everything above N defined as before, you can use the Matrix form referred to in the help page for NonlinearFit as follows, where again I use the default output from dsolve/numeric, i.e. procedurelist.

N1 := dsolve(eval(sys), {i(t), r(t), s(t)}, numeric, parameters = [beta, alpha]);
Vcompute := proc (t1, beta, alpha)
   N1('parameters' = [beta, alpha]);
   eval(i(t), N1(t1))
end proc;
TI:=Matrix([tV,iV]);
f:=proc(p, v, i) Vcompute(v[i,1],p[1],p[2])-v[i,2] end proc:
NonlinearFit(2,f,TI,initialvalues=);
#Notice that I have changed Vcompute slightly in order to make the implicit subtleties in Joe Riel's Vcompute explicit.
# t is the global t appearing in the output from e.g. N1(0.1). t1 is the formal parameter in the procedure. 

If the two animations are saved in p1 and p2, respectively, then

display(Array([p1,p2]));

will show them side by side and they will be played synchronously (assuming the same number of frames).

But since you already seem to have tried that, I wonder what you actually did.

The simplifications of 0^x to 0 and 0^0 t0 1 are done in Maple before any evaluations are taking place. This can be seen by enterering these two expressions surrounded by "unevaluation quotes":

'0^x';

resulting in 0 and

´0^0';

resulting in 1.

Compare with

'sin(0)';

and

sin(0);

There are a few other automatic simplifications made by Maple, e.g. 3+4 is simplified to 7.

For the system in question there are infinitely many critical points and worse: they are not isolated. The line given by x = -1, z = 0, and y in R\{0} consists of critical points.
(The image you provided clearly is of another system.)

restart;
f:=(x,y,z)->x*z/y-x-1;
g:=(x,y,z)->x*y-5*x*z+y;
h:=(x,y,z)->x*z^3/y+z^2/y+6*z;
sys:= diff(x(t),t) = f(x(t),y(t),z(t)),
      diff(y(t),t) = g(x(t),y(t),z(t)),
      diff(z(t),t) = h(x(t),y(t),z(t));
solve({f(x,y,z)=0,g(x,y,z)=0,h(x,y,z)=0},{x,y,z});
#That these are the only critical points can easily be verified (by hand or Maple) by solving the two first w.r.t. x,z:
solve({f(x,y,z)=0,g(x,y,z)=0},{x,z});
#The jacobian as a function can be found this way:
J:=unapply(VectorCalculus:-Jacobian([f(x,y,z),g(x,y,z),h(x,y,z)],[x,y,z]),x,y,z):
J(-1,y,0);
LinearAlgebra:-Eigenvalues(%);
#As expected at least one of the eigenvalues is zero. But since there is a positive eigenvalue you can state that all critical points are unstable.

In the help page for fsolve  and fsolve,details it says:

"For a general equation or system of equations, the fsolve command computes a single real root."

And under the explanation of the optional 'complex':

"complex
Find one root (or all roots, for polynomials) over the complex floating-point numbers."

You may try

fsolve(EQ1,T=300+100*I,complex);
fsolve(EQ1,T=300+100*I);

inspired by the results from

RootFinding:-Analytic((rhs-lhs)(EQ1),T=-5-5*I..400+400*I);

If you want a formula for the solution and not just a numerical solution, you could do either one of the following two:

##Solve with a general condition x(0)=a:

sol:=dsolve({diff(x(t),t$2) + 3*diff(x(t),t) + 11*x(t) = (t^2)*sin(t),x(0)=a});
eval(sol,a=0);
eval(sol,a=1);
#Both at once:
eval~(sol,[a=0,a=1]);

#Alternatively define sol as a function of a:
sol:=a->dsolve({diff(x(t),t$2) + 3*diff(x(t),t) + 11*x(t) = (t^2)*sin(t),x(0)=a});
sol(0);
sol(1);
#Both at once
sol~([0,1]);

If you want to (or have to) solve numerically, you can follow the advice given by PatrickT, i.e. use the parameters option for dsolve/numeric, but then you would have to give an initial value for the derivative too.

restart;
m:=Matrix(2, 2, {(1, 1) = phi*(w1*exp(mu)/(exp(mu)+1)+(1-w1)*exp(mu+eta[2])/(exp(mu+eta[2])+1)), (1, 2) = phi^2*(w1*exp(mu+tau[3])/(exp(mu+tau[3])+1)+(1-w1)*exp(mu+tau[3]+eta[2])/(exp(mu+tau[3]+eta[2])+1))*(1-w1*exp(mu)/(exp(mu)+1)-(1-w1)*exp(mu+eta[2])/(exp(mu+eta[2])+1)), (2, 1) = 0, (2, 2) = phi*(w1*exp(mu+tau[3])/(exp(mu+tau[3])+1)+(1-w1)*exp(mu+tau[3]+eta[2])/(exp(mu+tau[3]+eta[2])+1))});
expand~(m);
subs(exp(mu)=1/s1,exp(tau[3])=1/s2,exp(eta[2])=1/s3,phi=s5,w1=s6,%);
simplify(%);

Will this do what you want:

dsol := dsolve({eq1, eq2, x(0) = 0, y(0) = 0, D(x)(0) = Vx,D(y)(0) = Vy}, numeric, output = listprocedure,parameters=[Vx,Vy]);
X,Y:=op(subs(dsol,[x(t),y(t)]));
Xp:=proc(t,V,theta)
if not type([t,V,theta],list(realcons)) then return 'procname(_passed)' end if;
dsol(parameters=[V*cos(theta),V*sin(theta)]);
X(t)
end proc;
#And similarly you can define Yp if desired.