## A function of only one variable:

# R code to maximize h(x) = (x+47)*sin(sqrt(abs(256+x+47) )) + 512*sin(sqrt(abs(512-x-47)))
# over [-512, 512]

my.fcn.simp <- function(x){
output <- (x+47)*sin(sqrt(abs(256+x+47) )) + 512*sin(sqrt(abs(512-x-47)))
return(output)
}

my.n <- 100000

my.trial.xvalues <- cbind(runif(my.n, min=-512, max=512))
my.trial.hvalues <- my.fcn.simp(my.trial.xvalues)

my.max<-max(my.trial.hvalues)

position.max<-(1:my.n)[my.trial.hvalues==my.max]

location.max<-my.trial.xvalues[position.max]

print(paste("The maximum is ", round(my.max,3),
 " which occurs at x= ", round(location.max,3) ))

## To plot this function:

my.x <- seq(-512,512,length=50000)

plot(my.x,my.fcn.simp(my.x),type='l')




## A function of more than one variable:

# R code to maximize h(x1, x2) = (x1^2 + 4*x2^2)*exp(1 - x1^2 - x2^2)
# over [-3,3]^2

my.fcn <- function(x){
output<-(x[1]^2 + 4*x[2]^2)*exp(1 - x[1]^2 - x[2]^2)
return(output)
}

### NOTE:  This example sets up 'my.fcn' as a function of a vector 'x' which is (x1,x2).
### If the original function h was expressed as a function of two variables x and y, such as:
### h(x, y) = (x^2 + 4*y^2)*exp(1 - x^2 - y^2)
### then you could still write 'my.fcn' as a function of the two-dimensional vector 'x' and have 
### x[1] in place of x and x[2] in place of y.

my.n <- 1000

my.trial.xvalues <- cbind(runif(my.n, min=-3, max=3), runif(my.n, min=-3,max=3))
my.trial.hvalues <- apply(my.trial.xvalues, 1, my.fcn)
my.max<-max(my.trial.hvalues)

position.max<-(1:my.n)[my.trial.hvalues==my.max]

location.max<-my.trial.xvalues[position.max,1:2]

print(paste("The maximum is ", round(my.max,3),
 " which occurs at x1= ", round(location.max[1],3),
 " and x2= ", round(location.max[2],3)))



####
## Plotting this function:
####

my.fcn2 <- function(x,y){
output<-(x^2 + 4*y^2)*exp(1 - x^2 - y^2)
return(output)
}

my.x <- seq(-3,3,length=50)
my.y <- my.x
my.z <- outer(my.x,my.y,my.fcn2)

# One look:
persp(my.x, my.y, my.z)

# A different look:
persp(my.x, my.y, my.z, theta=90)

# With detailed tick marks:
persp(my.x, my.y, my.z, theta=90, ticktype = "detailed")