I know that something is very wrong here, but I absolutely cannot figure out what it is.

x<-c(-10:10)
xi<-cbind(rep(1,21),x)
y<-ifelse(x<=0,0,1)
nclusters = 2
mus = sample(x,nclusters)
sds = rep(1,nclusters)
betas <- matrix(rep(c(1,1),nclusters),ncol=nclusters,nrow=2)
pcs = rep(1/nclusters,nclusters)
posterior = matrix(nrow=21,ncol=nclusters)
ll = 0
for(q in 1:5){
  for(i in 1:nclusters){
    pyxcm = dnorm(y,mean=xi%*%betas[,i],sd=sds[i])
    pxcm = dnorm(x,mean=mus[i],sd=sds[i])
    pc = pcs [i]
    num = pyxcm*pxcm*pc
    denom = 0
    for(j in 1:nclusters){
      denom = denom + (dnorm(y,mean=xi%*%betas[,j],sd=sds[j])*dnorm(x,mean=mus[j],sd=sds[j])*pcs [j])
    }
    posterior[,i] = num/denom
  }
  for(i in 1:nclusters){
    inside1 = betas[1,i]*(y-xi%*%betas[,i])
    inside2 = betas[2,i]*(y-xi%*%betas[,i])
    temp1 = sum(inside1*posterior[,i])/sum(posterior[,i])
    temp2 = sum(inside2*posterior[,i])/sum(posterior[,i])
    ll = temp1*temp2
  }
  
  for(i in 1:nclusters){
      pcs[i] = mean(posterior[,i])
  }
  for(i in 1:nclusters){
    numer = sum(x * posterior[,i])
    denom = sum(posterior[,i])
    mus[i] = numer/denom
  }
  for(i in 1:nclusters){
    sds[i] = mean((x-xi%*%betas[,i])^2)
  }
  for(i in 1:nclusters){
    a = c(sum(y*x*posterior[,i])/sum(posterior[,i]),sum(x*posterior[,i])/sum(posterior[,i]))
    x2 = sum(x^2 * posterior[,i])/sum(posterior[,i])
    b = rbind(c(mus[i],x2),c(1,mus[i]))
    betas[,i] = solve(b)%*%a
  }
  for(i in 1:nclusters){
    temp = sum((y-xi%*%betas[,i])^2*posterior[,i])/sum(posterior[,i])
    sds[i] = ifelse(temp>.00001,temp,.1)
  }
}

numer = 0
denom = 0
for(i in 1:nclusters){
  numer = numer + ((xi%*%betas[,i])*dnorm(x,mean=mus[i],sd=sds[i])*pcs[i])
  denom = denom + (dnorm(x,mean=mus[i],sd=sds[i])*pcs[i])
}
yhat = numer/denom
plot(x,yhat)