@bliengme It appears that the book you are reading has the wrong result (or that you copied it incorrectly).
The notation Q[inf] seems to suggest that that should be the limit of Q(t) as t-> infinity.
That interpretation agrees with the stated result if r > 0.
But the same Q[inf] appears in P: That could be an error in copying from the book (or a misprint in the book).
Consider this where I test if Q__inf is just a constant and is the same appearing in both places:
 

restart;
P := t->Q__inf*r/(2 + 2*cosh(-r*t + b)); # Q[inf]
res:=int(P(s), s = 0 .. t); 
res1:=convert(res,exp);
res0:=Q__inf/(1+exp(b-r*t));
# A simple test:
eval(res1-res0,{r=7,b=3,t=5,Q__inf=1});
evalf[30](%);
# Numerical integration:
res2:=Int(P(s), s = 0 .. t); 
evalf[30](eval(res2-res0,{r=7,b=3,t=5,Q__inf=1}));

The results show that the book (or your copying) is wrong.
## NOTE: The two res0 and res differ by a constant since their derivatives are equal:
 

res0:=Q__inf/(1+exp(b-r*t));
##
diff(res0,t)-P(t); 
convert(%,exp);  # 0

@torabi SYSPDE is a set having 20 members. So far so good.
But take a look at e.g. SYSPDE[20].
I get:
-9.585649527*10^13*xi[3, 3] + 1.171579387*10^14*xi[3, 4]
     + 1.437847430*10^14*xi[4, 3] - 1.757369082*10^14*xi[4, 4]
     - 2.955249835*10^16*kappa[3, 3] + 3.149063384*10^15*kappa[3, 4]
     + 1.235128694*10^16*kappa[4, 3] - 4.260288677*10^13*kappa[4, 4] + Matrix(
    6, 6, [[-0.5521334135*10^15, 0., 0., 0., 0., 0.], [0., -0.5521334135*10^15
     - 2532.446978*kappa[3, 3] + 1266.223489*kappa[3, 4]
     + 1266.223489*kappa[4, 3] - 633.1117444*kappa[4, 4],
    -2849.002850*kappa[3, 3] + 1424.501425*kappa[4, 3],
    -2849.002850*kappa[3, 4] + 1424.501425*kappa[4, 4], 1266.223489*kappa[3, 3]
     - 2532.446978*kappa[3, 4] - 633.1117444*kappa[4, 3]
     + 1266.223489*kappa[4, 4], 0.], [0.,
    -2849.002850*kappa[3, 3] + 1424.501425*kappa[3, 4],
    -0.5521334135*10^15 - 3205.128206*kappa[3, 3], -3205.128206*kappa[3, 4],
    1424.501425*kappa[3, 3] - 2849.002850*kappa[3, 4], 0.], [0.,
    -2849.002850*kappa[4, 3] + 1424.501425*kappa[4, 4],
    -3205.128206*kappa[4, 3], -0.5521334135*10^15 - 3205.128206*kappa[4, 4],
    1424.501425*kappa[4, 3] - 2849.002850*kappa[4, 4], 0.], [0.,
    1266.223489*kappa[3, 3] - 633.1117444*kappa[3, 4]
     - 2532.446978*kappa[4, 3] + 1266.223489*kappa[4, 4],
    1424.501425*kappa[3, 3] - 2849.002850*kappa[4, 3],
    1424.501425*kappa[3, 4] - 2849.002850*kappa[4, 4], -0.5521334135*10^15
     - 633.1117444*kappa[3, 3] + 1266.223489*kappa[3, 4]
     + 1266.223489*kappa[4, 3] - 2532.446978*kappa[4, 4], 0.],
    [0., 0., 0., 0., 0., -0.5521334135*10^15]]) + 1.022469284*10^15*w[4, 3]
     - 1.533703927*10^15*w[4, 4] - 5.410566626*10^15*w[3, 3]
     + 7.157284986*10^15*w[3, 4]

You will notice that this is a sum of algebraic terms and a Matrix.
Try:
 

for i to nops(SYSPDE[20]) do i = op(i,SYSPDE[20]) end do;

whattype(op(9,SYSPDE[20]));

There may be similar problems for other SYSPDE[i].

@Dark Energy To plot diff(a(t),t)/a(t) use that diff(a(t),t) is given as the right hand side of ode1:
 

## Continuing with ode1 as defined previously:
ode1;
res:=dsolve({ode1,a(0)=1},numeric,parameters=[c,w,z]);
res(parameters=[1,2,3]); # Setting parameters of choice
plots:-odeplot(res,[t,a(t)^2-1],-1..10); # The easy one
## For the next you have to evaluate the right hand side at the parameter values:
plots:-odeplot(res,[t,eval(rhs(ode1),{c=1,w=2,z=3})/a(t)],0..10,thickness=3,labels=[t,diff(a(t),t)/a(t)]);

Note about the somewhat broken line from ca. 6 to 10:
rhs(ode1) is zero here, but this is a numerical computation, so roundoff happens. Thus some of the t values could result in imaginary numbers because of the square root. They can't be plotted, thus the holes.
 

eval(rhs(ode1),{c=1,w=2,z=3});
eval(%,res(8)); # An example where I get : 2.76925989236618*10^(-8)*I

In the help page for pdsolve under Parameters we find about 'build':
"(optional) try to build an explicit expression for the indeterminate function, no matter what the generality of the solution found"
(my emphasis).

After just a glance at table 3 (thanks to Tom Leslie) I find that what it says is:
"Proposed approach (approximate, Maple)".
Another Maple entry is:
"Proposed approach (Maple)".
Obviously, by "proposed approach" is meant methods proposed by the authors. So the difference between the two mentioned must be found by reading the article, which I haven't done.
There is no approximate Maple as opposed to just plain Maple.

@tomleslie

I agree with your three points:

1. Put the restart statement in its own execution group
2. Don't use 2D input
3. Don't use document mode

@Carl Love Even in document mode and with 2d input I get -1 as I should.

@Jodelkoenig I just tried in Maple 2017.3 with restart separated as Carl describes. I had no problems.
Which Maple 2017 do you have?
Mine is:
interface(version);
`Standard Worksheet Interface, Maple 2017.3, Windows 10, September 27 2017 Build ID 1265877`
and

kernelopts(version);
`Maple 2017.3, X86 64 WINDOWS, Sep 27 2017, Build ID 1265877`

@awass There is no file attached.

@vv Nice, clear, and short.

I think that in the line
 

P := indets( f, name ) minus { x }:

'minus {x}' should be removed. My reason is that the line defining Q:
 

# Solve for unknowns in g.
       Q := indets( g, name ) minus P:

is intended to get the unknown parameters in g and thus shouldn't include x.
## Actually you can replace P and Q by:

Q:={seq(a[k],k=0..n/2),r};

The question: You are allowing for several solutions. Does it ever occur?
 

@Josci95 You wrote that the code was supplied by Maple support. From the use of foldl, and the keen use of unevaluation quotes in things like 'parameters' = [ ':-a' = a, ':-b' = b ] I would have guessed that myself.
All very well, but why use foldl instead of plain add?

The worst part though, is that the code appears to be terribly inefficient.
In the uploaded worksheet below I have compared their approach to my traditional approach otherwise sticking to their settings. I used 1000 time values from the time interval 0..5. In the original there were only 5 values.
The essential difference between their approach and mine is that I have only one procedure Q. That procedure uses listprocedure output from dsolve.

solN := dsolve( { de, ic }, 'numeric','abserr'=1e-15, 'maxfun'=0, 'parameters'=[a,b],output=listprocedure ):
Xp,Yp:=op(subs(solN,[x(t),y(t)]));
###
Q:=proc(a,b) local j,ssq;
  if not [a,b]::list(realcons) then return 'procname(_passed)' end if;
  try
    solN(parameters=[a,b]);
    ssq:=add( (Xp(T[j])-X[j])^2,j=1..N)+add( (Yp(T[j])-Y[j])^2,j=1..N);
  catch:
    ssq:=10^15
  end try;
  ssq
end proc:

MaplePrimes19-01-08_ode_fit_test.mw
 

@Josci95 You could ask if there is any way of knowing that x(t) depends on the product a*b only without solving exactly as done here. In general that could be difficult if not impossible. But in this case we can see it by eliminating y(t) from the system.
Since we are pretending that we don't have data for y(t) that isn't a bad idea.
 

restart;
de := diff( x(t), t ) = a * y(t), diff( y(t), t ) = b * x(t);
ic := x(0) = 1, y(0) = 0;
###
diff(de[1],t);
ode_x:=subs(de[2]*a,%); ## depends on a*b only.
ics_x:=x(0)=1,D(x)(0)=0:  #The latter follows from ic and de[1].

@Rariusz Since theta(t) and diff(theta(t),t) are dependent on the groups c1, c2, and c3 only, the values of J, b, K, R, and L cannot be determined uniquely from data involving only theta(t) and diff(theta(t),t).
You will need data on i(t) too.

Example of the nonuniqueness: You can set R and L to whatever you want then J,b, and K are determined. Keep in mind though that the derivation of the 3rd order ode assumed that J, K, and L were nonzero.
 

## c-results:
Cr:= [c1 = 0.07085248546343231, c2 = 1072.44657764141, c3 = 91326.7791941768];
## sbs from my code:
sbs:={J = (L^2*c2 - L*R*c1 + R^2)/(L^3*c3^2), K = (L^2*c2 - L*R*c1 + R^2)/(c3*L^2),
    b = (L^2*c2 - L*R*c1 + R^2)*(L*c1 - R)/(L^4*c3^2)};
## Now pick L<>0 and R as you please:
eval(sbs,{R=1,L=2,Cr[]});
## or
eval(sbs,{R=10,L=7,Cr[]});

If you want you can even make b=0 for any L<>0 if you pick R = L*c1:
 

eval[recurse](sbs,{R=L*c1,Cr[]});

@Rariusz You write:
" Column 1 - time, column 2 - input signal, column 3 - position of rotor, column 4 - speed of rotor ".

Just for the record:
My interpretation related to eq1, eq2 is:

"input signal" =   V(t), "position of rotor" = theta(t), "speed of rotor"  = diff(theta(t),t).