Can reproduce the table 1 values. For table 2, the skin friction corresponds to the second solution but the fig 5 plot corresponds to the first solution, so something in the paper seems wrong. But the secret is to supply some approximate solutions that look like the solution you want.

Download ode.mw

To successively apply functions, use function composition with @. So expand@conjugate@i_ does your i_, then  conjugate, then expand. Edit: function multiplication at #2 now works

composition2.mw

Since the 19.58480220 is recognizable as having units of V^(-1), Vt must be in Volts, and so you want Vt from 0 to 0.0257 V. So you can solve for t and then use a parametric plot.

[Worksheet not displaying right now]

Download Vt.mw

Well, you are a moderator (note the symbol below your avatar), "a selected group of long-time MaplePrimes users with a reputation for positive contributions." - see https://www.mapleprimes.com/help/moderation.

There are some good reasons to edit others' posts; removing copyrighted content is one I've used. Another is removing duplicate copies of posted worksheets.

There are guidelines on the cited page. One of them is:

"In general, try to avoid making spelling or grammar changes unless they are particularly egregious, but feel free to edit message formatting if it messes up the display of the page." Since I don't think @acer's grammar was egregious, perhaps you are in violation :)

But seriously, this policy does not seem to have led to much difficulty.

@dharr OK, I see now what you want - here is how I would do it with laplace. Not sure about with Elziki

laplace.mw

solve(invEqs, {theta[1], theta[2], theta[3]}, explicit) gives two solutions, though they are many pages long. Once you have theta[1], theta[2] and theta[3], you can find sin(theta[1]) or various other derived quantities.

Your 7x7 Matrix has only Rank 6, so it must have Determinant zero, no matter what the values of k, Bi, omega etc are. If this isn't what your expect, you will need to check the construction of the Matrix. I didn't see any obvious reason for this, like two rows or columns the same.

Maple does have this as a builtin command.

CurveFitting:-PolynomialInterpolation(lxi, Lu, xi)

returns 

The dsolve command solves DAE systems. See the help ?dsolve,dae. There are five DEA methods,  rkf45_dae, ck45_dae, rosenbrock_dae, or mebdfi, and if you specify stiff=true, you will get the rosenbrock method. There is a brief desciption and some options described on that help page and the ?dsolve,dae_extension help page.

Download Number_solutions.mw

Here is one way. You may have to optimize plot parameters for the other plots.

Download sim_plot.mw

pdsolve can give a general solution, which can probably be worked up with BCs and IC in some cases to give an analytical solution, though not immediately in my version for the BCs you gave (incompatible?) :

Download pdsolve.mw

Numerical solutions are also possible, see ?pdsolve,numeric. [Edit - only for time and two space variables.]

To make b^1 a name, select it with the mouse, and then use the right-click (ctrl-click on mac)  menu to choose 2-D math -> convert to -> atomic variable.

atomic.mw

You can set up a Vector of r values, say rvals, at which you want the solution, and then use output=rvals. The output is a 2x1 Matrix, from which you can extract the information you want.

dsolvetable.mw

The option axes=none removes the axes and labels.