LinearAlgebra[ColumnOperation]  command is not suitable to solve your problem.

Here is a short way without any packages:

A:=<1,-1,2,-1; 2,-2,3,-3; 1,1,1,0; 1,-1,4,3>;  # Construction of a matrix

b:=<-8,-20,-2,4>;  # Construction of a vector

A1:=<b | A[..,2..4]>;  # Replacement of the first column of a matrix

 Addition.  More 2 examples:

A2:=<A[..,1] | b | A[..,3..4]>;  # Replacement of the second column of a matrix

A4:=<A[..,1..3] | b >;  # Replacement of the last column of a matrix

solid of revolution for  (a) about x=-1:

restart;

with(plots): with(plottools):

ax := spacecurve([-1, t, 0], t = -1 .. 6, color = red, thickness = 2, linestyle = dash):  # Axis of rotation

s1 := animate(plot3d, [[(x+1)*cos(t)-1, x, -(x+1)*sin(t)], x = 0 .. 2, t = 0 .. alpha, style = surface, color = green], alpha = 0 .. 2*Pi, frames = 100):

s2 := animate(plot3d, [[3*cos(t)-1, y, -3*sin(t)], y = 2 .. 5, t = 0 .. alpha, style = surface, color = green], alpha = 0 .. 2*Pi, frames = 100):

s3 := animate(plot3d, [[(x+1)*cos(t)-1, x^2+1, -(x+1)*sin(t)], x = 0 .. 2, t = 0 .. alpha, style = surface, color = green], alpha = 0 .. 2*Pi, frames = 100):

s4 := animate(plot3d, [[cos(t)-1, y, -sin(t)], y = 0 .. 1, t = 0 .. alpha, style = surface, color = green], alpha = 0 .. 2*Pi, frames = 100):

display(ax, seq(s||i, i = 1 .. 4), scaling = constrained, lightmodel = light3, axes = normal, orientation = [60, 70], labels = [x, y, z]);

For  revolution  (b) about y=6, everything can be done similarly.

My procedure is not as elegant as Douglas B. Meade's one, but it does exactly what OP wanted. In addition, it works in older versions Maple in which there is no element-wise operation ~ :

Cat:=proc(S)

local S1;

uses ListTools;

S1:=[Categorize((x,y)->is(x[2]=y[2]), [op(S)])];

{map(t->`if`(nops(t)=1,op(t), [[seq(t[i,1],i=1..nops(t))],t[1,2]]), S1)[]};

end proc:

Example :

M10 := {[0,[1,4,2,1]],[1,[2,1,4,2]],[2,[3,6,4,3]],[3,[3,6,4,3]],[4,[4,2,1,4]],[5,[4,3,6,4]],[6,[4,3,6,4]],[7,[5,8,6,5]],[8,[5,8,6,5]],[9,[5,8,6,5]],[10,[6,4,3,6]]}:

Cat(M10);

{[0, [1, 4, 2, 1]], [1, [2, 1, 4, 2]], [4, [4, 2, 1, 4]], [10, [6, 4, 3, 6]], [[2, 3], [3, 6, 4, 3]], [[5, 6], [4, 3, 6, 4]], [[7, 8, 9], [5, 8, 6, 5]]}

TheKingDanQc

Unfortunately I did not understand what you wrote. Here are the other 2 ways.

The easiest way to write such a procedure - it's just  to use a ready formula and write it as a procedure-function:

det := A::Matrix -> A[1, 1]*A[2, 2]*A[3, 3] - A[1, 1]*A[2, 3]*A[3, 2] - A[1, 2]*A[2, 1]*A[3, 3] + A[1, 2]*A[2, 3]*A[3, 1] + A[1, 3]*A[2, 1]*A[3, 2] - A[1, 3]*A[2, 2]*A[3, 1] :

Example of use:

A := <a1, b1, c1; a2, b2, c2; a3, b3, c3>;

det(A); 

The second method uses the definition of the determinant of any 3x3 matrix as the sum of the summands of the form  

sign((j-i)*(k-i)*(k-j))*A[i, 1]*A[j, 2]*A[k, 3]

 taken over all permutations of  {i, j, k} :

det1 := proc (A::Matrix)

local S, Rows, i, j, k;

S := 0; Rows := {1, 2, 3};

for i to 3 do

for j in Rows minus {i} do

for k in Rows minus {i, j} do

S := S+sign((j-i)*(k-i)*(k-j))*A[i, 1]*A[j, 2]*A[k, 3]

od; od; od;

S;

end proc:

The procedure  det1  is easily generalized to the case of a determinant of order 4, if you add 1 more loop.

The solution easy to get by  fsolve  command if we make a little preparatory work. The system can be reduced to two equations with two unknowns, if we subtract the second equation from the first equation. By  plots[implicitplot]  we find the initial approximations for  T[i]  and  T[r] , and the initial value for  H  we find by solving the first equation:

restart:

T[l]:= 20:

C[p]:= 3779:

rho:= 1026:

nu:= 0.004/60:

P:= 1000:

lign1:= H = 40.09 + 14.774*ln(T[i]-T[l])+(2.28+0.338*ln(T[i]-T[l]))*T[l]:

lign2:= H = -nu*rho*C[p]*ln(1-(T[i]-T[r])/(T[i]-T[l])):

lign3:= P = nu*rho*C[p]*(T[i]-T[r]):

plots[implicitplot]([lign1-lign2, lign3], T[i]=20..50, T[r]=20..50, color=red, thickness=2, gridrefine=3);

fsolve({lign||(1..3)}, {T[i]=30, T[r]=25, H=solve(eval(lign1,T[i]=30))});

An example of using numerical solution:

restart;

dsys := {diff(x(t),t)=y(t),diff(y(t),t)=-x(t),x(0)=1,y(0)=0}:

Sol:=dsolve(dsys,numeric, output=listprocedure);

x:=rhs(Sol[2]);

y:=rhs(Sol[3]);

fsolve(sin(x(t))+x(t)*y(t)^3-0.3=0, t=0..2);  # The solution of a nonlinear equation in a specific range

plot(sin(x(t))+x(t)*y(t)^3-0.3, t=0..2, color=red, thickness=2);  # Visualization

Edited. 

You can exactly to solve your equation  Eq  and then to define the function  U(x,t):

restart;

Eq:=diff(u1(z),z)=9/4*(84*(-z)^(11/2)-8*z^7+540*(-z)^(5/2)+324*z^4-324*z)*u1(z)/((z^2+3*sqrt(-z))*(2*(-z)^(3/2)+3)*(8*z^6+66*(-z)^(9/2)-189*z^3+216*(-z)^(3/2)+81)):

dsolve({Eq, u1(-2)=1});  # Exact solution

assign(%);

u0(z):=z^2/3+sqrt(-z):

U:=unapply(eval(u0(z)+t*u1(z), z=x/t),x,t);

plot(U(x,1), x=-2..0, color=red, thickness=2);  # The plot  U(x,t)  for  t=1

1) You forgot to call  CurveFitting  package. Instead  PolynomialInterpolation  should be  CurveFitting[PolynomialInterpolation]

2) f=-z  is the equation of degree 5 with 2 parameters. Even if there were no parameters, such equations in the general case can not be solved symbolically in terms of its coefficients (in radicals). See  wiki  and  wiki . In the best case, you can solve it numerically (if you specify parameter values).

da := [[1,m],[2,m2],[3,1],[4,2],[5,3],[6,4]];

f := CurveFitting[PolynomialInterpolation](da, z);

solution := solve(f=-z, z);

fsolve(eval(f=-z, [m=1,m2=2]), complex); # Numerical solution for specific values of parameters

Suppose you are looking for the lowest four-digit prime number whose sum of digits is divisible by 17. Write the program as a procedure without parameters:

restart;

Check:=proc()

local i, j, k, l, b;

for i from 1 to 9 do

for j from 0 to 9 do

for k from 0 to 9 do

for l from 0 to 9 do

b:=1000*i+100*j+10*k+l;

if isprime(b) and irem(i+j+k+l, 17)=0 then return i, j, k, l  fi;

od: od: od: od:

end proc:

Check();

                                       1, 0, 9, 7

Addition.  Using a procedure also has the advantage that all of the variables, to which you give any values in the procedure's body by default  will be local, and outside the body of the procedure they are ordinary symbols:

P:=proc()

local i, j;

i:=1; j:=2;

end proc:

i;  j;

                                 i

                                 j

See help on  Student[LinearAlgebra][GaussianEliminationTutor]  command.

If I understand you're looking for minimax optimization. DirectSeach  package provides a slightly better solution:

DirectSearch[GlobalOptima](max(abs(pa^2 + 2*pa*po - 0.44), abs(pb^2 + 2*pb*po - 0.10), abs(po^2 - 0.42), abs(2*pa*pb - 0.04)), [pa = 0.279246445735439, pb = 0.0728480349005200, po = 0.648114970425358]);

0.000685054125967331, [pa =0 .27924642145279, pb = 0. 07284794181838, po =0 .64811497042536], 925]

because  0.000685054125967331<0.000685109649562879

Sys := {diff(x(t),t) = y(t)-2*(x(t))^2, diff(y(t),t) = -8*(x(t))^3+4*x(t)*y(t)-5*(x(t))^4-(x(t))^5};

DETools[phaseportrait](Sys, [x(t), y(t)], t=-1..1, [[x(0)=1, y(0)=0]]);

R:=solve(f, x, explicit):

`*`(seq(x-R[i], i=1..4));

I pretty often use it. The operator  ~  seems to appear beginning with Maple 15 or Maple 16. It replaces not only  map  operator, but also  zip  operator. Here is an example in which you want to equate each element of the first list to the square of the corresponding element of second list:

X:=[x,y,z]:  A:=[a,b,c]:

zip(`=`, X, map(t->t^2,A))  # First way

X=~A^~2;  # Second way

                                 [x = a^2, y = b^2, z = c^2]

                                 [x = a^2, y = b^2, z = c^2]

This operator provides a more compact (short) syntax. I do not know his other advantages.

This works in one session until  restart  command is not used.

local I;

I:=2;

I+1;

                                      I := 2

                                         3