f(x,y) = 100(y-x^2)^2 + (1-x)^2

using PyPlot

x = -2:0.1:2
y = -2:0.1:2
X = repmat(x',length(y),1)
Y = repmat(y,1,length(x))
Z = map(f,X,Y)

fig = figure(figsize=(8,8))
ax = fig[:add_subplot](1,1,1, projection = "3d")
ax[:plot_surface](X, Y, Z, rstride=2,edgecolors="k", cstride=1, cmap=ColorMap("jet"), alpha=0.8, linewidth=0.25)
xlabel("x")
ylabel("y")

x = -2:0.01:2
y = -2:0.01:2
X = repmat(x',length(y),1)
Y = repmat(y,1,length(x))
Z = map(f,X,Y)
levels=float(0.0:0.25:10.0)

fig = figure(figsize=(10,6))
ax = fig[:add_subplot](1,2,1)
cp = ax[:contour](X, Y, Z, linewidth=2.0, levels=levels)
ax[:clabel](cp, inline=1, fontsize=10)
xlabel("x")
ylabel("y")

x = 0:0.01:2
y = 0:0.01:2
X = repmat(x',length(y),1)
Y = repmat(y,1,length(x))
Z = map(f,X,Y)
levels=float(0.0:0.25:2.0)

fig = figure(figsize=(10,6))
ax = fig[:add_subplot](1,2,1)
cp = ax[:contour](X, Y, Z, linewidth=2.0, levels=levels)
ax[:clabel](cp, inline=1, fontsize=10)
xlabel("x")
ylabel("y")

x = .95:0.001:1.05
y = .95:0.001:1.05
X = repmat(x',length(y),1)
Y = repmat(y,1,length(x))
Z = map(f,X,Y)
levels=float(0.0:0.01:.1)
ax = fig[:add_subplot](1,2,2)
cp = ax[:contour](X, Y, Z, linewidth=2.0, levels=levels)
ax[:clabel](cp, inline=1, fontsize=10)
xlabel("x")
ylabel("y")
tight_layout()

using SymPy

@syms x y
f(x,y) = 100(y-x^2)^2 + (1-x)^2
[diff(f(x,y), x); diff(f(x,y), y)]


[diff(f(x,y), x, x)  diff(f(x,y), x, y); diff(f(x,y), y, x) diff(f(x,y), y, y)]

A = [0 0 .5 1; .5 0 0 0; .5 1 0 0; 0 0 .5 0]

using JuMP
using Ipopt

m = Model(solver = IpoptSolver())
@variable(m, r[1:4] >= 0.)

@objective(m, Min, dot(r-A*r, r-A*r))

status = JuMP.solve(m)

@show rval = getvalue(r)

norm(rval-A*rval, Inf)